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THE  CLASSIFICATION  OF  MULTI-EDGE  SHAPES 
USING  AN 

AUTOREGRESSIVE  MODEL 
AND  THE 

KARHUNEN-LOEVE  EXPANSION 


by 

Ruth  DeBeurs  Kennett 

University  of  New  Hampshire,  December,  1985 

In  this  thesis  a  pattern  recognition  system  capable  of  classifying 
two  dimensional  shapes  with  multiple  edges  was  developed.  The  problem 
of  multiple  edge  classification  was  treated  as  an  extension  of  the  single 
edge  problem.  For  each  edge,  a  feature  vector  was  formed  from  the  para¬ 
meters  of  an  autoregressive  model  of  a  time  series  representing  the  shape 
of  the  edge.  The  dimension  of  these  feature  vectors  was  further  re¬ 
duced  by  the  use  of  a  transformation  based  on  the  Karhunen-Loeve 
expansion.  A  minimum  distance  classification  rule  was  used  to  classify 
an  input  transformed  feature  vector  according  to  the  nearest  class  mean 
in  the  transformed  feature  space. 

Two  boundary  sampling  methods  as  well  as  two  versions  of  the 
Karhunen-Loeve  transformation  were  investigated.  An  illustrative 
numerical  example  and  the  description  of  the  system  tests  are  provided. 
Using  an  equal  angle  boundary  sampling  technique  and  the  pre-whitened 
Karhunen-Loeve  transformation,  an  industrial  shapes  test  showed  100% 
correct  classification  results  with  an  average  classification  time  of 
1.27  seconds.  The  complete  Fortran  listings  of  the  routines  written  for 
this  system  are  included  in  the  Appendix  at  the  back  of  this  work. 
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CHAPTER  1 


INTRODUCTION 

The  purpose  of  this  thesis  is  to  classify  objects  with  multiple 
edges  using  the  method  of  classification  based  on  a  Karhunen-Loeve  trans¬ 
formation  of  the  autoregressive  model  parameters  which  represent  the 
shapes  of  the  boundaries  detected  in  thresholded,  digital  images  of  the 
objects. 

The  human  ability  to  recognize  a  particular  object  can  be  broken 
down  into  a  series  of  logical  steps.  For  instance  in  order  to  recognize 
a  cat,  a  person  looks  at  the  cat  and  unconsciously  collects  the  unique 
discriminating  features  of  the  animal.  Then  a  comparing  process  occurs 
where  the  just  collected  information  is  compared  to  the  vast  storehouse 
of  labelled  features.  A  match  occurs  when  the  newly  collected  features 
agree  with  the  features  remembered  as  those  belonging  to  a  cat. 

An  automatic  pattern  recognition  system  uses  the  same  principles. 

The  block  diagram  of  a  simple  pattern  recognition  system  is  shown  below 
in  Figure  1.1. 


Figure  1.1.1  Block  Diagram  of  a  Pattern  Recognition  System 


A  sensing  device  such  as  a  TV  camera  passes  the  raw  data  to  a  storage 
device.  In  the  storage  area  the  two-dimensional  image  can  be  digitized, 
thresholded,  segmented  and  preprocessed  in  any  way  to  enhance  the  dis¬ 
criminating  features  of  the  object.  Once  the  image  has  been  satisfactor¬ 
ily  preprocessed,  the  physical  features  or  the  structural  relationships 
(or  both)  can  be  selected  and  put  into  array  (also  called  pattern  vector 
or  feature  vector)  form  for  computer  manipulation.  At  this  point  the 
pattern  vector  can  undergo  further  processing,  or  it  can  go  directly  to 
a  classifier  or  categorizor.  Generally  the  classifier  has  been  pre¬ 
viously  trained  with  examples  of  expected  pattern  vectors  called  training 
samples  and  their  corresponding  labels.  Pattern  classification  occurs 
when  the  input  data  is  assigned  to  one  of  a  finite  number  of  categories 
or  pattern  classes. 

In  this  thesis  a  pattern  recognition  system  which  is  capable  of 
classifying  a  broad  set  of  objects  is  developed.  The  feature  vectors  of 
the  objects  are  formed  from  the  parameters  of  a  particular  linear  differ¬ 
ence  equation  called  the  autoregressive  (AR)  equation.  The  AR  equation 
has  been  shown  to  be  an  effective  model  for  the  shape  of  an  object's 
boundary  [Dubois,  1 984 j .  To  further  reduce  storage  space  and  classifi¬ 
cation  time,  the  AR  parameters  are  then  projected  onto  the  Karhunen- 
Loeve  (K-L)  coordinate  axes.  By  making  use  of  the  special  properties  of 
the  K-L  expansion,  the  dimension  of  the  feature  vector  can  be  further 
reduced.  The  result  of  this  work  has  been  the  development  of  an  efficient 
object  recognition  system  that  requires  very  little  training  data 
storage  space  and  no  extra  dedicated  hardware. 

The  list  of  applications  for  an  automatic  pattern  recognition  system 


is  almost  endless.  At  present,  recognition  systems  are  being  used  in 
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industry  to  inspect  and  identify  parts.  As  a  part  moves  down  an 
assembly  line  the  vision  system  can  inspect  for  quality  and  can  also  make 
sure  that  the  right  part  is  moving  down  the  line  at  the  proper  time. 
Electronic  suppliers  use  vision  systems  to  ensure  proper  placement  of 
circuit  board  components  [Beavers  and  Hubach,  1984].  Industrial  robots 
are  currently  in  use  for  parts  inspection,  sorting,  and  material  handling. 
Two  specialized  bin-picking  robots  are  already  finding  their  way  into 
factories;  one  is  the  GE  Bin  Vision  robot  and  the  other  is  the  ORS  (Object 
Recognition  System)  system  called  i-bot  1  [Edson,  1984].  Much  research 
is  being  done  for  the  autonomous  robot  where  the  ability  to  recognize 
obstacles  is  necessary  to  prevent  the  robot  from  walking  into  things  or 
falling  off  precipices. 

Other  areas  for  the  application  of  a  pattern  recognition  system  are 
in  the  military  and  the  medical  fields.  A  recognition  system  such  as  the 
one  developed  in  this  work  could  be  adapted  for  use  in  the  classification 
of  objects  in  remotely  sensed  images.  Other  military  applications  for  a 
pattern  recognition  system  include  automatic  undersea  searching  and 
target  classification.  In  the  medical  field  this  work  could  be  imple¬ 
mented  to  classify  cells  or  chromosomes. 

The  vision  systems  in  use  today  generally  require  a  great  amount  of 
training  data  in  order  to  classify  an  object.  A  reduction  in  the  amount 
of  data  needed  to  characterize  an  object  would  allow  more  storage  space 
for  the  training  data  of  additional  objects.  This  would  allow  the  vision 
system  to  be  able  to  recognize  more  objects  and  less  time  would  be  needed 
for  retraining.  The  recognition  of  three  dimensional  shapes  may  require 
the  storage  of  many  views  of  the  same  object.  This  is  only  possible  if 
the  amount  of  data  to  describe  one  view  is  minimal.  The  development  of 


a  recognition  system  that  could  classify  an  object  using  a  few  parameters 
would  be  beneficial  to  any  system  that  needed  to  classify  a  large  number 
of  objects. 

The  subsequent  chapters  will  introduce  and  then  describe  in  detail 
the  theory  and  implementation  of  the  pattern  recognition  system  developed 
in  this  work.  Chapter  2  introduces  some  feature  selection  techniques 
found  in  the  literature.  The  rest  of  the  chapter  then  focuses  on  the 
particular  autoregressive  model  and  the  parameter  estimation  algorithm 
used  to  produce  a  feature  vector.  Chapter  3  describes  the  classifi¬ 
cation  process  and  how  that  process  is  made  simpler  and  quicker  by  trans¬ 
forming  the  feature  vector  onto  a  set  of  Karhunen-Loeve  coordinate  axes. 
Chapter  4  discusses  the  actual  software  implementation  of  the  pattern 
recognition  theory  in  order  to  produce  a  working  system.  Chapter  5 
provides  a  numerical  example  of  the  system  theories  and  also  covers  the 
tests  of  the  system  and  their  results.  Chapter  6  concludes  this  work 
with  some  comments  on  the  test  results  and  suggestions  for  further 
improvements. 
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FEATURE  SELECTION  THEORY 


According  to  George  S.  Sebestyen  "Pattern  recognition  is  a  process 
of  decision  making  in  which  a  new  input  is  recognized  as  a  member  of  a 
given  class  by  comparison  of  its  attributes  with  the  already  known 
pattern  of  common  attributes  of  members  of  that  class"  [Sebestyen,  1962]. 
It  is  the  job  of  the  feature  selector  to  form  the  set  of  attributes  which 
are  most  representative  of  the  common  features  of  a  class. 

As  a  way  of  introducing  the  problems  encountered  in  the  feature 
selection  process,  some  algorithms  found  in  the  literature  will  be 
described  and  compared  in  the  first  section.  The  second  section  will 
present  the  feature  selector  used  in  this  work.  Included  in  that  section 
will  be  the  description  of  some  boundary  representation  schemes,  and 
also  the  theory  and  parameter  estimation  of  the  model  used  to  produce 


the  feature  vector. 


2.1  Feature  Selection  Background 


Feature  selection  algorithms  are  divided  into  three  basic  groups. 

The  first  group  is  comprised  of  heuristic  algorithms  which  are  based  on 
ad  hoc  rules  dependent  on  the  particular  objects  to  be  classified.  Since 
we  are  not  dealing  with  any  one  particular  type  of  object,  we  are  more 
interested  in  general  algorithms  which  can  be  applied  to  any  type  of 
object.  These  algorithms  form  the  following  two  categories:  the 
syntactical  and  the  mathematical.  The  main  difference  between  the  two 


types  of  algorithms  is  that  the  syntactical  algorithms  deal  with  the 
structural  information  of  the  image.  This  differs  from  the  mathematical 
approach  which  fits  an  analytic  model  to  the  physical  features  of  an 
image.  These  two  types  of  algorithms  are  described  in  further  detail  in 
the  following  paragraphs. 

2.1.1  Syntactical  Approach  to  Feature  Selection 

Just  as  the  meaning  of  a  word  is  often  context  dependent  so  is  the 
shape  of  an  object's  feature.  Thus  it  can  be  assumed  that  the  handle  of 
a  cup  will  be  curved  and  will  probably  join  the  cup  at  two  contact  points. 
The  syntactic  or  linguistic  approach  decomposes  complex  patterns  re¬ 
cursively  into  simpler  subpatterns  in  the  same  way  that  a  sentence  can  be 
decomposed  into  letters.  Patterns  are  described  by  their  basic  elements 
or  subpatterns  along  with  a  set  of  syntactic  rules  or  pattern  grammar. 

In  the  classification  of  patterns  described  by  a  syntactic  algorithm,  the 
classifier  performs  a  syntax  analysis  while  parsing  the  pattern  and 
answers  the  question  of  whether  or  not  the  pattern  belongs  to  the  language 
generated  by  the  grammar.  As  can  be  deduced,  this  method  is  very  complex 
to  implement  due  to  the  very  broad  range  of  objects  it  tries  to  classify. 
Much  research  is  being  done  in  this  area. 

2.1.2  Mathematical  Approach  to  Feature  Selection 

One  of  the  earliest  mathematical  feature  selection  techniques  is  the 
method  of  moments.  In  this  method  a  pattern  is  represented  by  its  two 
dimensional  moments  calculated  from  a  density  distribution  function  with 
respect  to  a  pair  of  axes  fixed  in  the  visual  field  [Hu,  1962].  These 


moments  and  other  moments  found  in  a  similar  manner  can  be  formed  into 
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linear  combinations  called  "moment  invariants."  These  are  so-called 
since  they  are  invariant  under  a  number  of  similarity  transformations. 

The  major  disadvantage  of  this  method  is  that  although  the  first  few 
moments  convey  significant  information  for  simple  objects,  they  fail  to 
do  so  for  more  complicated  ones.  Furthermore  the  computational  re¬ 
quirements  are  substantial  [Pavlidis,  1978]. 

Another  mathematical  feature  selection  technique  is  based  on  the 
values  of  the  Fourier  Transform  of  some  characteristic  of  the  pattern. 

This  technique  is  the  basis  of  algorithms  which  involve  the  calculation 
of  the  Discrete  Fourier  Transform  (DFT)  of  an  object's  characterizing 
function  which  is  commonly  the  boundary.  The  boundary  samples  can  be 
expressed  in  terms  of  tangent  angle  vs  arc  length.  They  can  also  be 
expressed  as  terms  of  the  complex  function  formed  from  the  boundary 
sample  position  where  the  x  axis  of  the  reference  plane  corresponds  to 
the  real  part  of  the  sample  and  the  y  axis  denotes  the  imaginary  part. 

The  main  disadvantage  of  the  resultant  Fourier  coefficient  shape  de¬ 
scriptors  is  that  not  all  of  the  coefficients  are  invariant  to  translation 
and  rotation  [.Granlund,  1972]. 

In  the  next  section  we  will  describe  the  theory  of  a  different 
mathematical  technique  for  feature  selection  based  on  the  autoregressive 
equation.  In  that  section  we  will  describe  how  the  AR  model  parameters 
can  be  estimated  and  how  they  are  formed  into  a  feature  vector. 

2.2  The  Autoregressive  Model 


The  autoregressive  model  is  a  probability  model  or  stochastic  model 
of  an  observed  time  series.  It  has  been  applied  in  the  areas  of  spectral 


analysis,  speech  recognition  and  transmission,  and  economic  forecasting. 
In  this  section  we  will  show  how  the  autoregressive  (AR)  model  can  be 
used  in  shape  recognition. 

2.2.1  Boundary  Representation  Schemes 

It  is  our  goal  to  compress  as  many  of  the  discriminating  features  of 
an  object  boundary(s)  into  as  few  AR  parameters  as  possible.  This  goal 
emphasizes  the  need  to  accurately  sample  the  boundary  to  catch  the  dis¬ 
criminating  features.  However  we  also  need  the  resultant  real  valued 
time  series  to  maintain  the  same  form  despite  changes  in  the  object's 
orientation,  size,  and  initial  sample  point.  If  the  object's  centroid 
is  known,  there  arise  many  possible  methods  of  sampling  the  boundary  to 
form  the  desired  time  series.  These  methods  will  be  described  below. 

2. 2. 1.1  Equi -Angle  Sampling.  Consider  that  the  object  centroid  is 
the  origin  of  the  cartesian  system.  Then  N  radius  vectors  can  be  pro¬ 
jected  from  the  origin  to  the  boundary  of  the  object  as  shown  below. 


Figure  2.2.1  Example  of  Equal  Angle  Boundary  Sampling  Method 


Let  us  make  the  requirement  that  the  angle  between  consecutive  radius 
vectors  is  constant  and  equal  to  2-rr/N .  As  Figure  2.2.1  shows,  there  are 
two  possible  time  series  produced  from  this  type  of  representation.  The 
first  is  the  time  series  r^ ,  r r^,  ....  r^  formed  from  the  lengths  of 
the  radius  vectors.  The  second  time  series  is  d-j ,  dg,  d^,  ....  d^  formed 
from  the  distances  along  the  boundary  between  consecutive  radius  vector 
to  boundary  intersection  points.  Using  this  equal-angle  sampling  method, 
the  portions  of  the  boundary  closer  to  the  centroid  are  sampled  relatively 
more  frequently  than  the  portions  farther  away  from  the  centroid.  Also 
in  the  case  of  highly  curved  shapes,  a  large  segment  of  the  boundary  may 
not  be  sampled  at  all.  For  example  in  Figure  2.2.2  the  tips  of  the 
letter  S  aren't  sampled. 


Figure  2.2.2  Example  of  Boundary  Sampling  Method  Showing  How 
Sections  of  Boundary  Can  Be  Skipped 


This  makes  the  letter  look  cut  off  at  the  ends  and  the  time  series  re¬ 


presentation  of  this  will  differ  from  the  one  when  the  letter  is  rotated 
so  that  the  tip  is  sampled. 

A  way  to  eliminate  these  problems  is  to  sample  the  boundary  so  that 
one  sample  is  always  the  same  distance  away  from  another  sample.  This 
method  is  described  in  the  next  section. 

2. 2. 1.2  Equal -Arc  Length  Sampling.  As  before,  consider  the  origin 
of  the  Cartesian  coordinate  system  to  be  the  centroid  of  the  object 
boundary.  Again  we  project  N  radius  vectors  from  the  centroid  to  the 
boundary.  However,  this  time  we  make  the  requirement  that  the  arc-length 
between  boundary  intersctions  remain  the  same  and  proportional  to  the 
total  length  of  the  boundary.  This  boundary  representation  is  shown 
below  in  Figure  2.2.2. 


Figure  2.2.3  Example  of  Equal  Arc  Length  Boundary  Sampling  Method 


As  can  be  seen,  the  boundary  is  now  uniformly  sampled.  As  before,  we 
can  form  two  time  series  from  this  method  of  boundary  sampling.  The 


i  d  rdjnutm1. 
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first  is  the  time  series  r^ ,  r ^  rg,  ....  r^  formed  from  the  radius 
vector  lengths.  The  second  is  the  series  0^,  Qg....,  eN  formed  from  the 
angles  between  r^ ,  r^. 

Out  of  the  four  methods,  two  were  actually  implemented.  Both 
implementations  formed  a  time  series  from  the  radius  vector  lengths.  One 
time  series  was  formed  using  the  equal  angle  sampling  method  and  the  other 
time  series  was  formed  using  the  equal  arc  length  sampling  method.  Let 
us  make  a  requirement  that  the  boundary  be  closed  and  that  it  may  not 
cross  over  itself.  We  then  have  a  discrete  time  series  that  is  stationary 
and  periodic. 

2.2.2  The  General  Autoregressive  Model 

In  the  general  AR  model  the  current  value  of  the  process  is  expressed 
as  a  finite,  linear  aggregate  of  previous  values  of  the  process  and  a  zero 
mean  white  noise  term  et  [Box  and  Jenkins,  1976].  If  we  let  the  values  of 

the  process  at  equally  spaced  time  intervals  t,  t-1 ,  t-2 . be  y^.,  yfc  1 

yt_2>  - then 

yt  =  Vt-i  +  Vt-2  +  •••  +  Vt-P  +  et  (2'2J) 

This  is  called  an  autoregressive  process  of  order  p.  In  Equation  2.2.1 
the  present  value  of  variable  y  is  regressed  on  the  previous  values  of 
itself.  Thus  the  model  is  autoregressive.  The  above  equation  can  be 
written  more  compactly  as 

*t  s  j,  +  et  <2-2-2> 

2 

The  model  contains  p+1  unknown  parameters  4>i ,  4>2»  ♦  ar,d  v  which  have 

2 

to  be  estimated  from  the  data.  The  additional  parameter  a  is  the 
variance  of  the  white  noise  process  e^. 


In  the  frequency  domain  the  AR  model  can  be  formed  into  a  transfer 
function  which  has  no  zeros.  Therefore  the  AR  model  is  also  called  the 


all -pole  model .  The  AR  model  is  also  linked  with  the  maximum  entropy 
method  (MEM)  in  the  area  of  spectral  analysis. 

2.2.3  The  Specific  Autoregressive  Model 

As  shown  by  Dubois  [Dubois,  1984]  an  effective  AR  model  for  use  in 
shape  recognition  is: 

m  _ 

r^  =  a  +  l  6jrt-j  +  ^wt  (2.2.3) 

J  1 

whose  form  was  originally  suggested  by  Kashyap  and  Chellappa  [Kashyap 
and  Chellapa,  1981].  In  this  form  we  have: 

rt  =  current  radius  vector  length 

rt_j  =  the  radius  vector  length  detected  j  time  intervals 

before  the  current  rt.  Collectively  these  are  called 
the  lag  terms, 
m  =  model  order 
/ffwj.  =  current  error 

{w^}  =  a  sequence  of  random,  independent,  zero  mean  samples 
(white  noise  sequence)  which  has  unit  variance,  i.e. 
E(w.u).)  =  1  if  i=j  and  =  0  otherwise 

'  J 

-3 -j , . . . , 9m  =  the  AR  parameters  or  lag  coefficients  to  be 
estimated  from  the  time  series 
(a, 6}  =  unknown  constants  to  be  estimated. 

Looking  back  at  the  general  autoregressive  model  of  Equation  2.2.1 
we  see  that  the  specific  model  used  for  shape  recognition  contains  two 
additional  terms  i  and  /e.  Comparing  the  two  equations,  /3wt  corresponds 
to  the  ef  term  of  the  general  model. 
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2.2.3. 1  Estimating  a  and  8  of  the  Specific  AR  Model.  Let  us  assume 
that  we  already  know  the  9.,  i  s  1,  ....  n.  Then  we  still  need  to  find 
estimates  of  a  and  3  in  Equation  2.2.3.  Now  we  will  define  the  average 
radius  vector  length  as  r  so  that 


N 


I 

t=l 


(2.2.4) 


where  N  is  the  number  of  radius  vectors  in  the  object  boundary  representa¬ 
tion.  We  can  rewrite  the  general  form  of  Equation  2.2.2  to  model  the  set 
of  differences  r^-r  as: 


or 


m 

l  e.(rt  ,-r)  +  /bu>. 
j=l  J  *  J 


(2.2.5) 


m  m 

rt  =  E  Vt-i  *  Stot  +  r  -  l  V 

r  j=l  J  1  J  1  j-1  J 

Substituting  in  Equation  2.2.2  for  rt  and  cancelling  the  common  terms  we 

find  the  constant  a  can  be  estimated  by  the  equation: 

m 

a  =  r  (1  -  l  9,)  (2.2.6) 

j-1  J 

Note  that  a  is  directly  proportional  to  the  mean  radius  vector  length  and 

thus  is  an  indicator  of  a  shape's  size. 

Now  that  we  have  an  estimate  for  a  we  still  need  to  find  an  estimate 

of  6.  Rewriting  Equation  2.2.3  to  separate  the  term  we  have 

m 


=  rt  '  a 


j=i  j  ‘-j 


(2.2.7) 


Squaring  both  sides  and  taking  the  expectation  of  the  results  we  get 
SEUt2]  -  E[(rt  -  „  -  ]  yt./] 

J  * 

Since  the  sequence  of  u>t  has  unit  variance  we  call  6  the  residual  variance 

and  it  can  be  estimated  as: 

^  N  m 


3  =  N  l  (rt  -  a  -  l  0,r  .y 

IN  t=1  t  j=]  J  t  J 


(2.2.8) 


Looking  at  Equation  2.2.8  we  see  that  8  is  roughly  proportional  to  the 
average  of  the  squared  radius  vector  lengths.  We  may  conclude  then  that 
a  possible  size  invariant  shape  descriptor  would  be  in  the  form  a //b. 


2. 2. 3. 2  Estimating  the  Coefficients  of  the  Specific  AR  Model  It  has 
been  shown  in  the  thesis  of  Dubois  [Dubois,  1984]  that  the  AR  parameters 
dj,  ...,  9  ,  ol  can  be  found  by  using  the  method  of  least  squares.  In  that 
work  the  model  parameters  were  chosen  so  as  to  minimize  the  expected  value 
of  the  squared  error  term  3  in  Equation  2.2.8.  In  the  experimental  work 

of  Dubois  the  equations  were  solved  using  matrix  methods  which  required 

3  2  2 

(m+1)  /3  plus  on  the  order  of  (m+1 )  operations  and  (m+1)  storage 

locations  [Makhoul ,  1975]. 

Due  to  the  special  properties  of  the  AR  equations  in  terms  of  the 
auto-correlation  function  the  parameters  of  the  model  can  be  estimated 
recursively. 

We  use  the  general  form  of  the  AR  equation  which  models  the  set  of 

differences  r^-r  in  Equation  2.2.5.  With  the  a  term  as  defined  in 

Equation  2.2.6  we  find  that  the  general  AR  model  of  Equation  2.2.2  based 

on  these  differences,  is  equivalent  to  the  specific  AR  model  of  Equation 

2.2.3  suggested  by  Kashyap  and  Chel lappa. 

If  we  multiply  both  sides  of  Equation  2.2.5  by  the  term  (r^  ^-r)  we 

get  [Box  and  Jenkins,  1976]: 

m 


(rt.k-r)(rt-?)  =  + 


Finding  the  expectation  of  both  sides  leads  to: 

m 


EC(rt_k-r)(rt_r-r)]  =  EC.Lei(rt.k-?)(rt.i-?)]  +  EC(  Vk"^^ 


The  term  E[(rt_k-r)/i3tot]  vanishes  for  k>0  since  the  sequence  is  un¬ 
correlated  with  the  (rt_k-r)  term.  For  k=0  this  term  is 

m  £  2  2 

•  Vpto.  +  Su*  ]  which  is  equal  to  8E[w.  ]  and  is  defined  as  o  . 
jtt  t  m 

Let  us  define 

Rk-j  =  EC<rt-k-f»rt-j-?>3  1  fr  J1(rt-k-?Hrt-j-?) 

which  is  the  covariance  function.  Note  that  R_^  =  R..  Rewriting  Equation 
2.2.9  and  changing  the  order  of  summations  on  the  right  side  we  get: 

"k  =  9lRk-l  *  92Rk-2  *  +  9mRk-m  for  k>0  <2-2-'°> 

We  now  have  the  autoregressive  difference  equation  in  terms  of  the  auto¬ 
correlation  function  of  the  sampled  data. 


E[  l  Yt- 

j=l  J  ' 


Equation  2.2.10  applies  for  any  model  order  of  the  AR  process.  Since 
there  is  a  different  set  of  parameters  for  each  model  order  we  will  denote 
the  kth  coefficient  in  an  AR  process  of  order  m  as  emk.  The  last  co¬ 
efficient  at  k  =  m  is  thus  8,.  This  term  is  known  as  the  partial  auto- 

mm 

correlation  coefficient  and  also  is  referred  to  as  the  mth  reflection  co¬ 
efficient.  We  can  now  form  a  set  of  equations,  one  for  each  model  order. 

k  ml  k-1  m2  k-2  m,n-l  k-m+1  mm  k-m 

for  k  =  1 ,  . .. ,  m  (2.2.11 ) 

These  are  called  the  Yule-Walker  equations.  In  matrix  form  these  equations 
may  be  written 


*o 

R1 

R2  ‘  1 

'  Rm-T 

"ml 

V 

R1 

R 

0 

R1  '  ' 

Rm-2 

em2 

R2 

R2 

R1 

1 

Ro  '  ' 

l 

•  R  , 
m-3 

1 

"m3 

» 

R3 

l 

1 

R  , 
m-1 

1 

Rm  9 
m-2 

l 

Rm  7  '  1 
m-3 

1 

1 

'  Ro  ^ 

’ 

i 

mm 

Rm_ 

The  above  equations  form  an  mxm  autocorrelation  matrix  which  is  Toeplitz 
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since  all  the  elements  along  the  diagonal  are  identical.  Also  the  matrix 
is  positive  definite  which  follows  from  the  positive  definite  property 
of  the  autocorrelation  function. 

If  all  the  terms  in  the  mxm  matrix  of  Equation  2.2.12  are  normalized 
to  R  ,  the  correlation  lag  at  t  =  0,  the  Yule-Walker  equations  can  be 
written  in  terms  of  the  resultant  normalized  autocorrelation  coefficients 
[Makhoul,  1975].  Thus  we  get 


T 

C1 

p  III 

l2 

Cm-1 

0ml 

C1 

C1 

1 

r  ill 

h 

c  ? 
m-2 

0m2 

C2 

C2 

. 

1 

C1 

1 

1 

1  '  1  ' 

i 

i 

Cm-3 

, 

l 

6  o 

m3 

1 

1 

C3 

1 

1 

5n-l 

Cm-2 

i  i  i 

Cm-3 

1 

l 

6 

mm 

1 

Cm 

(2.2.13) 


Ri 

where  C.  =  p-  i  =  0,. . . ,m 

o 

By  recognizing  that  the  right  hand  column  vector  contains  mostly  the  same 
elements  as  in  the  autocorrelation  matrix,  the  equations  can  be  solved 
using  simple  algebraic  reduction  techniques  as  shown  in  Box  and  Jenkins 
[Box  and  Jenkins,  1976]. 

The  solution  of  the  equations  reduces  to  the  following  recurrence 
equations : 


) 

mm 


m-1 

C  -  j  ■)  ,  .C  • 
m  ^  n»-l  *J  m-j 

nPi 

1  "  Vl  ,jCj 


j=l ,  . . . ,  m-1 


(2.2.14a) 


mj 


Vl ,j 


"mm  Jm-1  ,m-j 


(2.2.  I4t>; 


These  equations  are  attributed  to  Durbin  [Durbin,  I960].  They  show  how 
the  parameters  at  the  desired  model  order  are  calculated  using  the  para¬ 
meters  of  all  the  previous  model  orders.  To  find  the  mth  order  process 


17 


parameter  set,  the  parameter  sets  {e^},  ( e21  *  922  ^ » *  -  - »  *eml’9m2' 


0^}  are  all  calculated.  This  recursive  method  requires  only  2m  storage 


locations  and  irr  plus  on  the  order  of  m  operations  [Makhoul,  1975].  We 


can  now  augment  Equation  2.2.12  which  represent  the  k=l,...,m  equations 

m  o 


of  2.2.9  w'th  the  k=0  equation:  R 


=  E  9iRj  + 
j=l  J  J 


m 


The  result  becomes: 


Ro  R1  R2 


R1  Ro  R1 


R2  R1  Ro 


m 

‘Vi 
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R  Rm  ,Rm  ,  '  '  ' 

m  m-l  m-2 
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am 

91 

0 

02 
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1 
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1 

9m 

1 

0 

To  solve  the  above  equations  we  use  a  set  of  recurrence  relations  which 
produce  the  same  e^s  as  the  Durbin  algorithm,  except  for  a  sign  change 


on  the  parameters  due  to  assuming  the  general  AR  model  is  l  e.yt  .  =  e.. 

j=0  J  z 

The  algorithm  is  initialized  by: 


=  -R,/R0 


°i2  *  <>V>Ro 


(2.2.15a) 

(2.2.15b) 


the  recursion  for  k  =  2,3,...,m  is  given  by 


k-1 


°kk  =  '^Rk  +  .1 


k  +  Vl,jRk-jJ/o  k-1 


9ki  "  J k-1 ,  i  +  e kk° k-1 , k- i 


(2.2.15c) 
(2.2. 15d) 


This  form  of  the  recurrence  equations  has  two  advantages  over  the  original 

equations  2.2.14.  First,  the  correlation  coefficients  need  not  be 

2 

normalized  before  being  passed  to  the  algorithm.  Also  the  ok  term  can 
be  used  in  the  model  order  selection  process.  The  correct  model  order  is 
chosen  when  the  AR  model  has  arrived  at  the  best  fit  to  the  observed 
time  series. 

If  the  model  order  is  chosen  too  high  or  too  low,  the  AR  model  will 

not  provide  an  accurate  representation  of  the  sampled  process.  In  the 

2  2 
ak  equation  (Equation  2.2.15c)  we  see  that  ak  is  dependent  solely  on 

the  present  and  past  values  of  the  partial  correlation  coefficient.  The 
term  is  thus  an  indicator  of  the  error  and  has  the  property  that  it  de¬ 
creases  (or  remains  the  same)  as  the  order  of  the  model  increases 

[Makhoul,  1975].  A  suitable  model  order  for  the  process  can  be  chosen 

2  2 

when  the  change  in  the  error  term  i.e.,  ak  -  a  k_-j  is  below  a  specified 
tolerance  level. 


2.2.4  The  Autoregressive  Model  Parameters  as  Shape  Descriptors 

The  shape  descriptors  form  a  feature  vector  which  is  a  string  or 
set  of  numbers  which  represent  the  features  or  characteristics  of  the 
shape.  Since  the  physical  features  of  a  shape  are  invariant  to  rotation, 
translation,  and  scale,  so  should  the  numerical  shape  descriptors  be  in¬ 
variant  to  rotation,  translation  and  scale.  Also  the  numerical  shape 
descriptors  should  be  invariant  to  where  the  sampling  process  begins. 

From  the  previous  section  we  have  estimated  the  following  AR  parameters: 

i  -  l,...,m  estimated  using  the  recursive  algorithm  of  Equation 
2.2.15,  a  as  in  Equation  2.2.6  and  3  from  Equation  2.2.8.  These  parameters 
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in  vector  form  [0-j ,  Qg,  ...,^1  a//e]^  produce  a  feature  vector  which  has 
the  desired  invariant  properties  [Dubois,  1984].  However,  it  should  be 
noted  that,  in  using  the  equi -angle  boundary  sampling  method,  the  rotation 
and  starting  point  invariance  property  only  applies  when  the  change  in 
rotation  or  starting  point  is  a  multiple  of  2tt/N.  If  N  is  sufficiently 
large,  this  no  longer  is  a  problem. 


CHAPTER  3 


CLASSIFICATION  OF  THE  TRANSFORMED  FEATURE  VECTOR 

Using  the  theory  from  Chapter  2,  it  is  possible  to  form  a  feature 
vector  from  the  set  of  AR  parameters  which  describe  the  shape  of  an 
object's  boundary.  If  we  had  a  bin  full  of  the  exact  same  object  and  if 
we  found  the  feature  vector  of  the  outer  boundary  of  each  object  using 
the  same  method,  we  would  find  that  most  feature  vectors  would  be  slightly 
different  from  the  rest.  This  is  due  not  only  to  minor  variations  in 
rotation,  i.e.  not  exact  multiples  of  2it/N,  but  also  to  pixel  quantization 
error  and  round-off  errors  in  the  system.  Considered  as  a  class  however, 
all  these  feature  vectors  will  be  roughly  the  same.  It  is  then  possible 
to  use  a  classification  scheme  based  on  the  feature  vector  of  the  samples 
of  each  class.  Some  of  these  classification  schemes  as  found  in  the 
literature  will  be  described  in  the  next  section.  It  will  become  apparent 
that  there  is  usually  a  need  for  further  processing  of  the  feature  vector 
in  order  to  correctly  classify  the  shape  of  the  object's  boundary.  This 
is  especially  true  for  objects  with  very  similar  shape.  Thus  the 
Karhunen-Loeve  ( K-L )  expansion  will  be  introduced  and  two  implementations 
of  the  K-L  expansion  will  be  described.  The  last  section  of  this  chapter 
will  be  devoted  to  the  particular  classification  rule  used  in  the 
experimental  work  of  this  thesis. 

3.1  Background  of  Classification  Theory 

i  . 


The  pattern  classification  problem  is  essentially  the  problem  of 
partitioning  the  feature  space  so  as  to  assign  each  possible  feature 


vector  or  point  in  that  space  to  the  proper  pattern  class.  The  classi¬ 
fication  techniques  in  the  literature  form  three  groups  depending  on  how 
much  a  priori  information  about  the  feature  vector  to  be  classified  is 
known.  The  first  group  comprises  the  deterministic  techniques  wherein  no 
assumptions  are  made  on  the  statistical  properties  of  a  pattern  class. 

When  the  statistical  properties  of  the  feature  vector  are  known,  then  it 
is  possible  to  use  the  group  of  statistical  classification  techniques.  A 
third  group  which  allows  incomplete  knowledge  of  the  statistical  properties 
of  the  feature  vectors  is  the  group  of  trainable  classifiers.  These  three 
methods  of  classification  will  be  described  in  further  detail  in  the 
next  pages. 


3.1.1  Deterministic  Classification  Techniques 

This  group  of  classification  techniques  depends  on  the  formulation  of 
mathematical,  deterministic,  discriminant  functions  which  are  used  to 
partition  the  feature  space.  These  functions  assume  that  the  feature  vector 
components  are  deterministic  quantities.  The  hyperplane  classification 
method,  the  third  technique  used  in  the  thesis  of  Dubois  [Dubois,  1984], 
is  an  implementation  of  a  deterministic  discriminant  function. 

The  general  form  of  a  discriminant  function  is  as  follows  [Fu,  1968]. 
Let  u»i  ,^2*  •  •  •  be  designated  as  the  c  possible  pattern  classes  to  be  re- 
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be  the  feature  vector  where  y.  represents  the  ith  feature.  The  dis¬ 
criminant  function  D.(Y)  associated  with  pattern  class  w.,  j=1,...,  c 

J  J 

is  such  that  if  the  feature  vector  Y^  is  in  class  then  the  value  of 
D-(Y)  must  be  the  largest  (or  possibly  the  smallest)  of  all  the  D,(Y)  for 

*  J 

j?h'.  The  discriminant  functions  form  decision  boundaries  between  regions 
associated  with  class  u.  and  u.  which  can  be  expressed  as: 

*  J 

Di(Y)  -  Dj(Y)  =  0 

An  example  of  the  discriminant  function  for  a  two  class  problem  is  shown 
in  Figure  3.1.1. 


Figure  3.1.1  Example  of  a  Discriminant  Function  for  a  2  Class  Problem 


An  important  class  of  linear  deterministic  classifiers  uses  as  the 
classification  criterion  the  distance  between  the  unlabelled  feature 
point  and  the  labelled  points  of  the  feature  space.  A  discriminant 
function  for  such  a  so  called  minimum  distance  classifier  is: 


»Tr*v 


where  N  is  the  number  of  samples  of  the  jth  class  and  R.  .  is  the  ith 

J  >  n , 1 

component  of  the  nth  sample  of  the  jth  class.  Equation  3.1.1  can  be 

implemented  more  simply  without  loss  of  information  as: 

o  N  d  ? 

MY)2  =  l  l  (Y-j  ”  Rj  n  i)  J  ■  1 . c 

J  n=l  i=l  1  J*0’1 

Classification  of  the  unlabelled  feature  vector  Y  is  based  on  the  smallest 


class  distance  D.{Y).  The  performance  of  the  minimum  distance  classifier 

J 


is  dependent  on  whether  or  not  the  training  sample  points  form  tight 
clusters  that  are  well  separated  from  each  other.  In  classification 
terminology  we  wish  for  the  training  samples  to  form  clusters  such  that 
their  intraset,  or  within-class,  distances  are  small  and  their  interset, 
or  between-class,  distances  are  large. 


3.1.2  Statistical  Classification  Techniques 


In  the  deterministic  classification  techniques,  it  is  assumed  that 
the  feature  vector  components  are  deterministic  quantities.  However,  in 
many  cases  the  noise  effect  due  to  large  variations  in  the  feature  vector 
components  cannot  be  neglected.  If  we  assume  that  for  each  class  the 
probability  density  function  of  the  feature  vector  is  known  and  the  pro¬ 
bability  of  the  ith  class  occurring  is  known,  then  the  classification 
problem  becomes  a  statistical  problem.  The  task  of  the  classifier  is 
then  to  minimize  the  probability  of  misclassification  [Fu,  1968]. 

The  statistical  approach  to  the  pattern  recognition  problem  is 
generally  based  on  the  Bayes  rule  which  presents  an  optimum  measure  of 
classification  performance  [Tou  and  Gonzalez.  1974].  The  classifier  that 
uses  the  Bayes  rule  is  called  a  Bayes  classifier.  One  basic  Bayes 
classifier  is  the  implementation  of  the  following  discriminant  function: 


24 


¥ 


fi- 


v 


°j(Y)  =  P(Y/^J)p(oJj.)  j  = 


A  pattern  vector  Y_  is  assigned  to  class  u.  if  for  that  class 

J 

D.(Y)  >  D-(Y)  for  all  i ^ j .  The  main  problem  with  the  statistical 

J  ' 

classification  techniques  is  that  of  estimating  the  probability  density 
function  characterizing  each  of  the  different  pattern  classes.  This 
implies  that  the  computational  requirements  of  these  techniques  will 
probably  be  greater  than  for  the  deterministic  classification  techniques. 


3.1.3  Trainable  Pattern  Classifiers 


The  trainable  pattern  classifiers  produce  discriminant  functions  by 
means  of  iterative  learning  algorithms.  These  classifers  can  use  either 
the  deterministic  or  statistical  techniques  or  a  combination  of  both  to 
produce  the  desired  discriminant  functions.  The  main  use  of  a  trainable 
classifier  is  when  the  required  information  for  optimal  pattern  classi¬ 
fication  is  only  partially  known.  During  its  operation  the  classifier 
learns  the  needed  information.  If  the  learned  information  gradually 
approaches  the  true  information,  then  the  decisions  based  on  the  learned 
information  will  eventuallly  approach  the  optimal  decision  as  if  all 
the  information  required  is  known  [Fu,  1968].  The  trainable  classifiers 
are  divided  into  two  groups:  the  supervised  classifiers  and  the  un¬ 
supervised  classifiers.  The  difference  between  the  two  groups  is  that 
the  supervised  classifiers  have  available  the  correct  classification  of 
the  observed  patterns  while  the  unsupervised  classifiers  don't.  Again, 
the  computational  demands  of  the  trainable  classifier  algorithms  are 
generally  greater  in  comparison  to  the  deterministic  techniques. 
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3.2  The  Karhunen-Loeve  Expansion 

It  is  our  goal  to  classify  objects  as  quickly  and  accurately  as 
possible.  Thus  it  is  desired  that  the  actual  classifier  should  be  a 
simple  deterministic  minimum  distance  classifier.  As  was  stated  pre¬ 
viously,  for  the  minimum-distance  classifier  to  give  accurate  results, 
it  is  necessary  for  the  feature  space  to  form  tight  clusters  well 
separated  from  each  other.  Also  if  the  feature  vector  dimension  can  be 
reduced,  the  classification  time  will  decrease.  In  this  section  we 
will  introduce  the  Karhunen-Loeve  expansion  and  its  properties  which  allow 
us  to  achieve  quick,  accurate  classification  of  feature  vectors  of  re¬ 
duced  dimension.  Two  implementations  of  the  K-L  expansion  will  be  de¬ 
scribed  and  finally  the  actual  minimum-distance  classifier  used  in  this 
work  will  be  introduced. 

3.2.1  Basic  Theory  of  Karhunen-Loeve  Coordinate  Axes 

Let  us,  first  of  all,  introduce  the  Karhunen-Loeve  expansion.  We 
know  that  a  periodic  function  can  be  expanded  into  a  Fourier  series: 

oo 

x(t)  =  l  xn  exp(jnw0t)  (3.2.1) 

n=-°° 

where  is  the  angular  frequency  and  the  xn  are  the  Fourier  coefficients. 
A  similar  expansion  can  be  formed  for  a  non-periodic  function  x(t)  where 

oo 

x(t)  *  l  ynup(t)  0<t<T  (3.2.2) 

n=l 

Suppose  we  use  only  D  discrete  sampling  points  to  approximate  the  time 
function  during  the  time  interval  0<t<T,  then  x(t)  becomes  the  column 
vector  x.  Also  the  generalized  Karhunen-Loeve  expansion  of  Equation  3.2.2 
becomes 


where  the  vectors  u  form  a  deterministic  orthornormal  basis  set.  The 
— n 

basis  vector  is  D-dimensional  and  represented  as 


(3.2.4) 


The  derivation  of  the  Karhunen-Loeve  coordinate  axes  from  the  K-L 
expansion  is  described  next,  following  the  method  of  Devijver  and  Kittler 
[Devijver  and  Kittler,  1982].  We  start  by  making  an  approximation  to  x 
of  Equation  3.2.3.  Let's  call  this  approximation  x,  and  let  it  have 
dimension  d  where  d<D  i.e. 


* =  l  yn± 


(3.2.5) 


The  problem  is  now  of  finding  the  vectors  u^  such  that  the  mean  square 
error  (MSE)  between  x  and  its  approximation  x  is  minimized.  Let  this 


error  be  called  e  where 


£  =  E{  (x-x_r(x-x)} 


(3.2.6) 


and  the  expected  value  E{}  is  over  all  samples  in  a  class.  Substituting 


in  the  expansions  for  x  and  x  we  get 


'  EUnl,  y^)' 


(3.2.7) 


Due  to  the  orthonormality  property  of  the  basis  vectors  we  know  that 


^n  hti  '5nm 


(3.2.8) 


where  6  is  the  Kronecker  delta  and  equals  1  if  n=m  and  0  otherwise.  Th 
nm 

error  of  Equation  3.2.7  reduces  to 


Since  we  still  want  to  minimize  the  error  in  terms  of  the  vectors  u^, 
let  us  find  yp  in  terms  of  u^.  Since  n  is  a  dummy  variable  we  can 


rename  it  and  rewrite  Equation  3.2.3  as 

D 

x  =  l  y  u 
—  S  m-m 
m- 1 

We  can  now  multiply  both  sides  of  Equation  3.2.10  by  u^  to  get 


t  t 

li  x  =  u  Y  yu 
— n  —  — n  ■'m-m 

m- 1 


(3.2.10) 


(3.2.11) 


Using  the  orthonormality  property  of  Equation  3.2.8  we  find  that  the 
right  side  of  Equation  3.2.11  will  be  zero  except  when  m=n.  Thus  we  get: 


yn  =  ^n  — 


n  =  1 ,. . .  ,D 


which  can  be  substituted  back  into  Equation  3.2.9  i.e. 

D  t  t 

£  I  V*  2i  V 

n=d+1  ^  ^ 


(3.2.12) 


(3.2.13) 


Since  the  u^  basis  vectors  are  deterministic  the  order  of  summation  and 


expectation  can  be  interchanged  to  get: 

0  t  t 

e  =  l  VtE{:*  * 

n=d+l  ^  ^ 

We  can  denote  the  expectation  term  as  w  where 

=  E{x  xS 

then  the  MSE  of  the  approximation  becomes 


u  t 

l  “nV 


n=d+l 


^  T^n 


(3.2.14) 


(3.2.15) 


(3.2.16) 


Since  we  wish  to  minimize  the  error  e  at  the  same  time  as  maintaining  the 
constraint  in  Equation  3.2.8,  we  need  to  find  a  form  of  the  error  equation 
where  the  desired  constraint  is  explicit.  To  find  this  we  use  the  method 
of  Lagrange  multipliers.  First  we  can  rewrite  the  constraint  as: 


u  “  u  -  1  =  0 
— n  — n 


(3.2.17) 


Multiplying  both  sides  of  Equation  3.2.17  by  a  constant  An  doesn't  change 
the  equality.  Thus  Equation  3.2.16  can  be  written  as  a  function  of  u^  as: 


D  t  D 

9(0  =  l  JL, V-  -  l  \  ‘  ^ 

^  n=d+l  ^  ^  n=d+l  n  ^  ^ 


(3.2.18) 


If  the  above  equation  is  differentiated  with  respect  to  u^  and  the  result 


is  set  equal  to  zero  we  get  the  following  condition: 

t”-An1)Hn  =  0 


(3.2.19) 


Thus  the  optimal  vectors  u^  are  the  eigenvectors  of  the  matrix  v  and  the 
constants  An  are  the  corresponding  eigenvalues.  Since  v  is  a  DxD  symmetric 
matrix,  it  will  have  at  most  D  linearly  independent,  orthogonal  eigen¬ 
vectors  and  D  real  eigenvalues  [Lanczos,  1956]. 

We  still  need  to  find  the  d  eigenvectors  out  of  the  set  of  D  in  order 
to  form  the  approximation  x  in  Equation  3.2.5.  From  Equation  3.2.19  we 
know  that  ^  =  A  u^.  Substituting  this  into  Equation  3.2.16  and  using 
the  orthonormal ity  condition  again,  we  get  the  error  in  terms  of  the 


eigenvalues  of  • p ,  i .e. 


D  .  D 

l  !LV  =  l  a. 

u  — n  — n  L  > 


n=d+l 


n=d+l 


(3.2.20) 


Here  we  can  see  that  to  minimize  e,  the  sum  on  the  right  should  include 
the  smallest  eigenvalues  of  the  matrix  .  This  implies  that  x  should  be 
formed  from  the  eigenvectors  corresponding  to  the  d  largest  eigenvalues 
of  >  i .e. : 


x  '  y" 

Where  the  eigenvalues  \n  corresponding  to  satisfy 
'1  -  x2  -  '3  -  •  ••  -  Ad  -  ■  ”  -  ‘D 


(3.2.20) 


(3.2.21) 


The  vectors  u„  form  what  are  called  the  Karhunen-Loeve  coordinate  axes. 
— n 
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It  has  just  been  shown  that  in  using  the  eigenvectors  corresponding 
to  the  largest  eigenvalues  of  the  matrix  v,  the  mean  square  error  be¬ 
tween  the  original  vector  and  its  approximation  is  minimized.  Also  it 
can  be  shown  [Devijver  and  Kittler,  1982]  that  the  coefficients  yn  of 
the  expansion  are  uncorrelated.  This  suggests  how  the  K-L  coordinate 
axes  can  be  used  to  reduce  the  dimension  of  the  input  vector  by  de- 
correlating  the  components  and  removing  the  axes  which  do  not  convey  much 
information  [Kittler,  1975]. 


mixture  set  containing  elements  of  many  classes  ^ ,  i  =  1 . c  where 

c  is  the  number  of  pattern  classes.  If  we  assume  that  the  samples  of 
each  class  are  normally  distributed  about  the  class  mean  vector,  then 
each  class  of  samples  can  be  completely  described  by  its  mean  vector,  m^ 
and  its  covariance  matrix  C.  [Tou  and  Gonzalez,  1974]  where  C-  is  de¬ 


fined  as: 


ci  =  E  C  ( xrHIi )  ( ) t  U 


i  =  1 ,. . .  ,c 


ni 


*  zr-  l  [(x-mi)(x-nv.)t]  i=l,...,c 

ni  k=l  1  1 

where  is  the  number  of  samples  of  the  ith  class  and 

=  ElX-j  i  =  1 ,. . .  ,c 


(3.2.24) 


i  "l 


i  =  1 , . . .  ,c 


(3.2.25) 


Using  this  additional  information  we  can  form  a  within  class  scatter  matri 
S  [Devijver  and  Kittler,  1982]  which  is  the  average  of  all  the  class  co- 

W 


variance  matrices  i.e.: 


s  =  y  p,c. 


i=l 


(3.2.26) 


where  P.  is  the  apriori  probability  of  a  sample  of  the  mixture  set  being 
from  the  ith  class.  If  a  sample  is  equally  likely  to  come  from  one  class 
as  any  other  then: 

(3.2.27) 


P.  =  1 
1  c 


and  the  matrix  S  becomes: 
w 


n . 


(3.2.28) 


sw  =  f  njJ,  <*fk  -  EfXiik  - 
Hence  the  columns  of  the  matrix  U  in  Equation  (3.2.22)  can  be  formed  from 


the  eigenvectors  of  Sw  ordered  by  the  decreasing  value  of  their 
corresponding  eigenvalues. 
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It  can  be  shown  that  the  population  entropy  can  be  minimized  by 
forming  the  matrix  V  from  the  eigenvectors  of  matrix  Sw  corresponding  to 
the  smallest  eigenvalues  [Devijver  and  Kittler,  1982].  A  good  measure 
of  intraset  dispersion  is  the  entropy  function  given  by: 

H  =  -E[log  p(y_)] 

where  p(y)  represents  the  class  probability  density  function  [Tou  and 
Gonzalez,  1974]. 

It  is  noted  that  the  within  class  scatter  matrix  Sw  is  symmetric  and 
positive  semi -definite,  thus  the  eigenvalues  will  be  real  and  greater  than 
or  equal  to  zero.  The  eigenvectors  point  in  the  direction  of  the 
principle  axes  of  a  hyperellipsoid  whose  shape  is  defined  by  the  co- 
variance  matrix  and  the  corresponding  eigenvalues  determine  the  length  of 
these  axes.  In  other  words,  the  eigenvalues  represent  the  variance  of  the 
data  in  the  direction  of  the  corresponding  eigenvectors  in  the  feature 
space. 


3.3  The  Optimal  Karhunen-Loeve  Coordinate  System 

There  are  many  variations  to  the  Karhunen-Loeve  coordinate  system  in 
the  literature.  These  variations  are  specially  tailored  to  their  particular 
type  of  pattern  recognition  problem  and  essentially  differ  in  the  choice  of 
the  K-L  system  generating  matrix.  These  techniques  are  inferior  to  the 
next  method  to  be  described  in  the  sense  that  they  do  not  utilize  the 
class  mean  information  in  an  optimal  manner  [Kittler,  1975].  In  this 
section  the  theory  of  the  optimal  Karhonen-Loeve  coordinate  axes  will  be 
described. 


Jy. 
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Generally  the  most  important  discriminating  information  is  contained 
in  the  class  mean  vectors  m^ .  Let  us  then  introduce  a  new  matrix  based 
on  the  class  mean  variances  referred  to  as  the  between  class  scatter  matrix 
Sfa  [Devijver  and  Kittler,  1982]  where: 


sb  =  l  P.(m--m)(m.-m)t 

D  i=]  1  -1  -1 


(3.3.1) 


Here  m^  is  the  class  mean  vector  as  defined  in  Equation  3.2.25  and  m  is 
the  average  of  the  class  mean  vectors  i.e.: 

1  c 

m  =  Em.  (3.3.2) 

c  i=1  1 

When  there  are  c  classes,  the  space  spanned  by  the  class  mean  vectors  will 
be  at  most  c-1  dimensional  [Devijver  and  Kittler,  1982].  Here  we  see  that 
by  using  matrix  Sb  as  the  K-L  coordinate  system  generating  matrix  we  may 
achieve  a  more  efficient  way  of  reducing  the  dimension  of  the  pattern 
vector . 

The  algorithm  as  suggested  by  Kittler  in  [Devijver  and  Kittler,  1982] 
is  described  next.  We  first  set  out  to  decorrelate  the  noise  on  the 
pattern  vectors.  From  the  previous  section  we  found  that  this  is  done  by 
the  transformation  matrix  V  in  Equation  3.2.22  formed  from  the  eigen¬ 
vectors  of  S.  We  then  would  like  to  normalize  the  variances  of  the 
w 

noise  components  of  the  pattern  vectors  to  unity.  To  do  this  we  form  a 
matrix  a  whose  diagonal  elements  are  the  eigenvalues  corresponding  to  the 
eigenvectors  of  V.  Combining  the  two  matrices  we  achieve  the  desired 
results  by  using  the  transformation  matrix  B  where 

B  =  V  A"1/2  (3.3.3) 


(3.3.4) 


Similarly  we  can  decorrelate  and  prewhiten  the  between  class  scatter 
matrix  Sb  of  Equation  3.3.1  to  form  Sb'  i.e. 

Sb'  =  BtSbB  (3.3.5) 

We  find  the  optimal  coordinate  system  by  finding  the  eigenvalues  and  eigen¬ 
vectors  of  the  matrix  Sb'.  The  optimal  feature  extracting  transformation 
matrix  is  now  W  where 

W  =  UA"1/2V  (3.3.6) 

where  V  is  the  matrix  whose  columns  are  formed  from  the  eigenvectors  of 
$b‘  in  order  of  the  descending  values  of  the  corresponding  eigenvalues. 
Finally  the  transformed,  reduced  dimension  feature  vector  of  Equation 
3.2.22  becomes 

£  =  W*x  (3.3.7) 

where  x^  is  the  original  m+1  dimensional  feature  vector. 

3.4  The  Classifier 


By  using  either  Equation  3.2.22  or  Equation  3.3.7  we  can  reduce  the 
dimension  of  the  feature  vector  and  also  decorrelate  their  components. 
Thus  the  resultant  transformed  sample  space  should  have  fairly  well  de¬ 
fined  boundaries  between  the  classes.  For  classification,  all  we  need  is 
a  simple  minimum  distance  classifier  as  described  in  Section  3.1.1.  The 
classifier  implemented  in  this  work  is: 


where 


l  (bk  -  yik) 

k*l 


2 


D.j  =  total  distance  for  the  ith  pattern  class 
bk  =  the  kth  dimension  of  the  transformed  unlabelled  (input) 
feature  vector 

yik  =  average  kth  component  for  the  ith  class  of  transformed 


training  vectors. 
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The  classification  decision  is  chosen  in  favor  of  the  smallest  class 
2 

distance  .  This  means  that  the  input  vector  is  categorized  as  be¬ 
longing  to  the  pattern  class  whose  average  sample  it  is  nearest.  Once 
the  classifier  has  been  trained  the  only  information  it  needs  to  classify 
an  input  object  is  the  transformation  matrix  and  the  mean  sample  trans¬ 
formed  feature  vector  of  each  expected  object  boundary. 

In  this  chapter  we  have  shown  how  the  dimension  of  the  feature  vector 
can  be  reduced  by  projecting  the  vector  onto  a  reduced  system  of  K-L 
coordinate  axes.  It  was  also  shown  how  the  class  mean  information  can  be 
optimally  used  to  maximize  the  interset  distances.  The  next  chapter  will 
present  the  necessary  hardware  and  software  algorithms  to  implement  the 
theory  of  this  and  the  previous  chapter. 


i'ttl 
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CHAPTER  4 


CREATING  THE  PATTERN  RECOGNITION  SYSTEM 


This  chapter  deals  with  the  realization  of  the  operational  pattern 
recognition  system  and  the  details  of  the  algorithms  used  in  this  system. 
The  first  section  of  the  chapter  will  describe  the  system  hardware  and 
how  the  original  scene  is  digitized  for  computer  use.  The  second 
section  discusses  the  thresholding  technique  used  on  the  digital  re¬ 
presentation  of  the  scene.  In  the  third  section  the  algorithms  used  to 
extract  the  raw  data  from  which  the  features  can  be  selected,  are  de¬ 
scribed.  The  fourth  section  describes  how  to  extract  and  use  the 
additional  information  provided  by  the  multiple  edges  of  an  object.  The 
last  section  describes* and  provides  the  flow  charts  of  the  software 
written  for  this  pattern  recognition  system. 


4.1  Object  Sensinc 


We  have  been  assuming  that  the  machine  that  will  do  the  necessary 
data  manipulation  for  feature  selection  and  classification  is  the  digital 
computer.  However  the  computer  doesn't  have  eyes  to  see  the  object!  We 
now  describe  how,  by  the  use  of  the  proper  hardware  and  special  software, 
the  computer  can  receive  data  just  as  if  it  could  see. 


4.1.1  Description  of  the  System  Hardware 


The  main  hardware  elements  of  the  object  recognition  system  developed 
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RS-170  standard,  2)  the  Image  Technology  Inc.  (ITI)  IP-512  series  image 
acquisition  board  containing  the  AP-512  Analog  Processor  and  the  FB-512 
Frame  Buffer,  3)  the  VAX  11-730  digital  computer  and  4)  the  Panasonic 
TR-920M  monitor  (black  and  white).  See  Figure  4.1.1  for  the  system 
design. 

The  TV  camera  scans  the  scene  and  produces  an  analog,  video  signal 
which  reproduces  the  variations  in  the  scene  intensity.  This  analog  signal 
is  sent  to  the  ITI  image  acquisition  board  where  it  is  flash  converted 
by  the  circuits  of  the  AP-512  Analog  Processor  into  a  data  stream  of  8  bit 
bytes.  Each  byte  is  a  pointer  to  one  of  256  elements  of  a  programmable 
look-up-table  (LUT).  If  the  byte  represents  the  number  8  then  the 
contents  of  the  8th  element  of  the  LUT  is  put  into  an  appropriate  location 
of  a  256- K  byte  RAM  called  the  FB-512  Frame  Buffer.  The  portion  of  the 
frame  buffer  memory  actually  used  in  the  system  is  called  the  active 
video  window.  This  window  is  a  480  row  by  512  column  two  dimensional 
array  of  8  bit  numbers  (called  grey  levels)  which  stores  the  digitized 
representation  of  the  complete  scene  as  viewed  by  the  TV  camera.  Each 
element  of  this  digitized  representation  is  called  a  picture  element  or 
pixel .  Since  the  TV  camera  starts  its  scanning  cycle  at  the  top  leftmost 
corner  of  the  scene,  the  origin  of  the  active  video  window  is  also  at  the 


top  leftmost  corner  (as  seen  on  the  monitor).  Figure  4.1.2  shows  the 
origin  of  the  video  screen  coordinate  system  of  the  active  video  window. 
The  x  values  increase  to  the  right  and  the  y  values  increase  downward. 
The  contents  of  the  frame  buffer  are  flash  converted  from  digital  to 
analog  for  real-time  viewing  on  the  monitor.  By  the  use  of  the  ITI 


software  applications  package,  the  Fortran  programs  running  on  the  VAX 


Pattern  Recognition  System  Hardware 


Figure  4.1.2  Active  Video  Window  Coordinate  System 
computer  can  access  the  contents  of  the  frame  buffer  and  the  hardware 


functions  of  the  analog  processor. 

4.1.2  The  Digital  Image 

In  order  to  provide  the  computer  with  data  in  the  form  that  it  can 
accept,  the  scene  containing  the  objects  is  sampled.  What  we  have 
stored  in  the  frame  buffer  then  is  a  digital  representation  of  the 
original  continuous  scene.  Each  pixel  of  this  representation  is  a  rec¬ 
tangular  element  which  has  an  8  bit  value  which  quantizes  the  image 
intensities  to  any  one  of  256  gray  levels. 

As  a  direct  result  of  this  quantized  representation  formed  from 
these  rectangular  pixels  there  arise  two  problems  which  must  be  con¬ 
sidered  in  the  implementation  of  the  system.  First,  the  shape  of  the  pixel 
causes  the  distance  measurements  (radius  vector  lengths),  based  on  pixel 
location,  to  vary  with  changes  in  object  rotation  or  translation.  In  the 
digital  representation  4  increments  in  the  x  direction  have  the  same 
physical  distance  as  5  increments  in  the  y  direction.  This  is  referred  to 
as  the  pixel  aspect  ratio. 
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To  eliminate  the  problem  due  to  the  rectangular  pixel  shape,  we 
introduce  a  pixel  aspect  correcting  factor  of  5/4  to  the  x  distance  values 


in  our  calculations. 

A  second  problem  that  arises  is  the  effect  of  the  quantization  pro¬ 
cess  on  the  original  scene.  Instead  of  the  smooth  edges  which  exist  in  the 
original  scene  the  representati ve  edges  are  uneven  and  sometimes  jagged. 
Also,  the  number  of  pixels  on  the  object  boundary  will  vary  with  slight 
changes  in  boundary  position.  Thus  we  see  variations  in  the  representa¬ 
tion  of  similar  objects  if  they  differ  slightly  in  two  dimensional 
translation  or  rotation.  These  variations  are  collectively  referred  to 
as  quantization  noise.  The  significance  of  this  noise  can  be  reduced  by 
decreasing  the  size  of  the  pixels  or  by  increasing  the  relative  size  of 
the  objects  in  the  representation. 

4.2  Preprocessing  of  the  Digital  Image 

Once  we  have  acquired  the  digital  representation  of  a  scene,  we 
generally  like  to  preprocess  the  representation  so  as  to  make  the 
feature  selection  algorithms  simpler  and  quicker.  The  bag  of  tricks 
which  allows  us  to  perform  this  preprocessing  comes  under  the  heading  of 
Digital  Image  Processing.  However,  due  to  the  time  constraint  on  a 
pattern  classifier,  the  preprocessing  of  the  image  is  kept  to  an 
essential  minimum.  In  this  section  we  will  describe  the  thresholding 
done  to  the  digital  image  in  order  to  simplify  the  boundary  detection 
algorithms. 


4.2.1  Image  Thresholdin 


In  order  to  implement  a  quick  and  simple  boundary  detection  scheme, 
the  boundaries  of  the  objects  must  be  distinct  from  the  background  and 
they  must  be  uniform  in  intensity.  A  simple  way  to  produce  these  desired 
properties  is  to  transform  the  image  into  a  binary  picture  which  contains 
objects  of  uniform  intensity  against  a  background  of  a  different 
intensity.  The  process  of  producing  such  a  binary  image  is  called 
thresholding. 

To  perform  this  thresholding  process  we  use  a  single  programmable 
look-up-table  on  the  ITI  image  acquisition  board.  On  the  board  there  are 
actually  2  groups  of  4  input  LUTs,  and  3  groups  of  4  output  LUTs.  The  2 
input  groups  allow  two  possible  camera  inputs  and  the  3  output  groups 
allow  for  the  three  color  groups  (RGB  of  a  color  system).  Each  of  the 
four  LUTs  in  any  group  can  be  preprogrammed  or  pre-loaded  before  actual 
use.  Then  during  execution,  the  program  can  switch  LUTs  and  instantly 
transform  the  scene.  For  our  purposes,  however,  we  don't  need  to  trans¬ 
form  the  scene  from  moment  to  moment,  but  we  do  need  to  threshold  the 
image  before  it  is  stored  in  the  frame  buffer.  To  do  this  the  256 
locations  of  the  chosen  input  LUT  are  loaded  with  one  of  two  grey- 
levels.  The  locations  below  a  certain  location  (called  the  threshold) 
are  loaded  with  one  grey  level,  for  example  0,  and  the  locations  on  and 
above  the  threshold  are  loaded  with  a  different  grey  level,  for  instance 
255.  Imagine  that  we  have  a  light  object  on  a  dark  background  and  a  grey 
scale  of  0  (black)  to  255  (white).  If  the  LUT  has  been  preloaded  as 
described,  then  as  the  scene  is  initially  viewed,  the  data  bytes  are 
immediately  converted  to  either  one  of  the  two  designated  grey  levels. 


The  light  areas  of  the  scene  where  the  intensity  is  above  the  threshold 
are  stored  as  grey  level  255  in  the  frame  buffer  and  the  dark  areas 
below  the  threshold  are  stored  as  grey  level  0. 
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4.3  Feature  Selection 

We  now  assume  that  we  have  a  binary,  thresholded  representation  of 
the  original  scene  stored  in  the  frame  buffer.  We  still  need,  however, 
to  form  the  time  series  from  the  radius  vector  lengths  of  a  sampled 
boundary  in  order  to  calculate  the  components  of  a  feature  vector. 

To  do  this,  an  object  must  first  be  detected.  Once  an  object  is 
detected,  the  boundaries  need  to  be  traced  and  the  boundary  centroids 
estimated.  When  we  have  a  centroid  pixel  position,  the  radius  vector 
lengths  can  be  calculated  and  the  corresponding  time  series  can  be 
produced.  It  is  only  when  the  complete  time  series  is  formed  that  the 
AR  parameters  can  be  calculated  using  the  algorithm  described  in  chapter 
two. 

In  the  next  pages  we  will  describe  the  various  techniques  used  to 
provide  the  raw  measurements  from  which  the  features  can  be  selected. 

4.3.1  Search  for  a  Boundary  Starting  Pixel 

This  pattern  recognition  system  allows  the  existence  of  more  than 
one  object  in  the  scene  viewed  by  the  camera.  Also  the  objects  may  have 
more  than  one  closed  edge.  In  order  to  recognize  these  multiple  edges, 
it  is  necessary  to  scan  the  entire  scene  stored  in  the  frame  buffer. 

To  speed  up  this  scanning  process  not  all  of  the  pixels  are  read. 
The  search  for  object  edges  begins  at  the  origin  of  the  active  video 


window  in  the  frame  buffer.  The  search  continues  along  the  constant  y 
row  where  only  the  x  values  which  are  multiples  of  a  previously  prompted 
for  increment  are  read.  Once  the  end  of  the  row  has  been  reached,  the 
search  continues  at  the  beginning  of  the  next  y  row  which  is  a  multiple  of 
the  increment.  For  example,  if  the  incremental  step  is  8,  the  following 
pixels  on  the  first  row  are  read:  (0,0),  (8,0),  (16,0),  ....  The  pixels 
which  are  read  on  the  next  row  are:  (0,8),  (8,8),  (16,8),  and  so  on. 

The  search  continues  in  the  same  manner  until  a  transition  occurs  in¬ 
dicating  an  object  has  been  detected  or  the  bottom  right  corner  of  the 
scene  has  been  reached. 

A  transition  is  detected  when  the  difference  between  the  present 
pixel  grey  level  and  the  previous  grey  level  is  equal  to  the  difference  in 
intensity  between  the  object  and  the  background  (in  our  case  255).  Due  to 
the  manner  of  hopping  across  the  rows  of  the  frame  buffer  in  search  of 
objects,  there  is  no  guarantee  that  the  pixel  where  the  transition  is  first 
encountered  is  an  edge  pixel.  Thus  the  search  backtracks  pixel  by  pixel 
on  the  same  row  until  the  transition  indicating  the  actual  edge  is 
detected.  This  pixel  where  the  boundary  is  first  detected  is  referred 
to  as  the  start  pixel. 

4.3.2  The  Turtle  Boundary  Detector 

When  an  object  in  its  two-dimensional  representation  is  of  uniform 
intensity  and  its  edge  doesn't  contain  spurious  gaps,  then  it  is  possible 
to  use  a  "turtle"  boundary  detection  algorithm.  In  this  work  a  modified 
version  of  the  turtle  algorithm  used  in  the  work  of  Dubois  [Dubois,  1984] 
is  implemented.  In  this  algorithm  the  boundary  between  the  object  and  its 
background  is  followed  in  a  clockwise  manner  by  a  symbolic  device  referred 


to  as  the  "turtle".  Once  the  boundary  start  pixel  has  been  found  as  de 
scribed  in  the  previous  pages  the  turtle  proceeds  to  trace  the  edge 
according  to  the  following  two  rules  [Duda  and  Hart,  1973]: 

1)  If  the  turtle  is  in  the  object  it  will  turn  left 
and  take  a  step. 

2)  If  the  turtle  is  in  the  background  it  will  turn 
right  and  then  take  a  step. 

The  turtle  stops  when  it  has  returned  back  to  the  start  pixel. 

Since  a  great  portion  of  the  classification  time  is  spent  turtling 
around  the  boundaries  of  the  objects  in  the  scene  it  is  advantageous  to 
speed  up  the  algorithm.  To  do  this,  the  turtle  algorithm  of  Dubois 
[Dubois,  1984J  was  modified  in  the  following  ways: 

1)  The  direction  of  the  turtle  movements  is  governed 
by  decisions  based  on  logical  flags  indicating  past 
and  present  turtle  position. 

2)  All  subroutine  calls  have  been  eliminated  in  pre¬ 
ference  for  a  quicker  straight  through  code. 

3)  As  each  boundary  pixel  is  detected,  the  x  and  y 
locations  are  temporarily  stored  and  counted  for 
later  use  by  the  radius  vector  calculation 
algorithms. 

4)  Also,  as  each  boundary  pixel  is  detected,  its  grey  level 
intensity  is  changed.  This  marking  of  the  boundary 

is  to  prevent  repeat  turtling  of  the  edge  as  the 
entire  frame  buffer  is  scanned.  Also  the  boundary 
marking  provides  identification  for  later  determination 
of  which  edges  are  inside  other  edges.  An  example  of 
the  scene  with  marked  edges  is  shown  in  Figure  4.3.1. 
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Figure  4.3.1  Example  of  a  scene  showing  marked  edges 
Note  that  the  edges  are  marked  with  the  intensity  level  which 
corresponds  to  their  order  of  detection.  For  instance  the  boundary 
closest  to  the  top  left  corner  is  detected  first  and  is  thus 
marked  with  the  grey  level  intensity  1. 


4.3.2. 1  The  Centroid  of  the  Boundary.  The  center  of  the  two- 
dimensional  representation  of  the  object  can  be  approximated  by  calculat 
ing  the  average  of  the  boundary  pixel  locations.  Thus  the  xc  and  yc 
coordinates  or  the  boundary  center  are: 
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where  N^  is  the  number  of  boundary  pixels  and  xn  and  yn  are  the  x  and  y 

coordinates  of  the  nth  boundary  pixel.  'The  advantage  of  finding  the 

boundary  centroid  is  that  despite  changes  in  the  object's  rotation  and 

translation  the  relative  position  of  the  centroid  remains  the  same.  The 
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centroid  position  can  then  be  used  as  the  origin  of  a  set  of  reference 
axes  needed  for  the  determination  of  the  boundary  intersection  radius 
vector  lengths. 

4.3.3  Formation  of  the  Time  Series  from  the  Radius  Vector  Lengths 

With  the  position  of  the  boundary  centroid  known,  the  two  types  of 
time  series  formed  from  the  radius  vector  lengths,  as  described  in 
chapter  two,  can  be  produced.  The  type  of  time  series  formed  depends  on 
the  relationship  between  each  radius  vector  of  the  boundary  representation 
Two  methods  based  on  two  types  of  relationships  are  explored  in  this 
work  and  described  below: 

4.3.3. 1  Equi -angle  Method.  In  this  method  the  intersections  of  the 
radius  vectors  with  the  boundary  are  such  that  the  angles  between  the 
radius  vectors  are  all  constant  and  equal  to  2iT/Nr  radians,  where  Nr  is 
the  number  of  radius  vectors.  The  algorithm  actually  implemented  is  a 
modified  version  of  the  algorithm  written  by  Dubois  [Dubois,  1984].  The 
main  modification  of  the  original  algorithm  is  the  use  of  the  boundary 
locations  stored  in  an  array  during  initial  boundary  detection.  This 
eliminates  the  additional  time  of  turtling  around  the  boundary  for  a 
second  time.  However  the  trade  off  is  the  additional  temporary  storage 
requirement  for  the  boundary  locations.  Another  modification  is  the  use 
of  straight  through  code  which  also  increases  the  speed  of  forming  the 
time  series. 

The  basic  steps  needed  to  form  the  desired  time  series  using  this 
method  are  described  next.  Initially  the  origin  of  a  cartesian  co¬ 
ordinate  system  is  shifted  to  the  centroid  of  the  object  boundary,  see 
Figure  4.3.2. 


Figure  4.3.2  Example  of  the  Equi-Angle  Radius  Vector 
Boundary  Intersection  Method 

The  radius  vectors  are  then  imagined  to  project,  equally  spaced  by  2ir/Nr 
radians,  from  the  origin  to  the  object  boundary.  Each  section  formed 
by  a  pair  of  radius  vectors  is  referred  to  as  a  sector.  Note  that  in 
determining  the  intersection  of  the  radius  vector  with  a  boundary  pixel, 
we  are  only  concerned  with  the  relative  position  of  the  pixel  with 
respect  to  the  radius  vectors  of  the  sector.  Thus  only  the  magnitudes 
of  the  slopes  of  the  first  quadrant  of  the  reference  plane  need  to  be 
calculated.  These  slopes  are  calculated  once  and  stored  in  an  array 
before  the  actual  search  for  boundary  intersections  begins.  The  slope 
of  the  jth  radius  vector  is  calculated  using 

slope(j)  *  5/4  *  tan(2Tr/Nr*j)  j=0 . (Nf/4)  -  1  (4.3.2) 

The  5/4  factor  is  introduced  in  Equation  4.3.2  to  deal  with  the  rectangul 
pixel  shape  and  to  prevent  the  factor  from  having  to  be  included  in  each 
pixel  slope  calculation  during  the  search. 


The  search  for  radius  vector  boundary  intersections  begins  with  the 
start  pixel.  The  slope  of  the  imaginary  radial  line  drawn  from  the  start 


pixel  (xs>  ys)  to  the  centroid  pixel  (xc,  yc)  is  given  as: 


slope  = 
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(4.3.3) 


This  slope  value  is  compared  with  the  slope  values  stored  in  the  array 
and  the  appropriate  sector  is  found.  The  two  slope  vectors  of  both  edges 
of  the  sector  form  reference  lines  for  calculating  the  relative  position 
of  the  pixel  in  the  sector.  After  the  start  pixel  has  found  its 
sector  and  the  relative  position  in  the  sector  has  been  noted,  the  second 

pixel  next  to  the  start  pixel  is  checked.  Its  relative  position  is 

noted  with  respect  to  the  reference  vectors  also.  If  the  distance 
measurements,  showing  the  relative  pixel  position,  change  signs  from 
one  pixel  to  the  next,  then  it  can  be  assumed  that  the  most  recently 
checked  pixel  is  either  on  or  it  has  crossed  a  radial  line.  The  radius 

vector  length  between  this  pixel  (x,y)  and  the  centroid  can  be  calculated, 

i  .e. 


rt(1)  -  /  (<x-xc)  *  5/4)2  +  (y-yc)2  (4.3.4) 

This  radius  vector  length  becomes  part  of  the  time  series  in  the  order  in 
which  it  is  detected.  This  process  repeats  till  all  the  boundary  pixels 
have  been  checked  and  the  start  pixel  has  been  encountered  once  again. 

The  details  of  the  implementation  of  this  method  are  described  in  the  work 
of  Dubois  [Dubois,  1984]. 


4. 3. 3. 2  Equi-Arc  Length  Method.  In  this  method  the  intersections 
of  the  radius  vectors  with  the  boundary  are  such  that  the  arc  length 
between  intersections  is  constant  and  equal  to  a  certain  number  of  pixels. 
In  the  implementation  of  this  method  the  constant  arc  length  (A), 
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measured  in  pixels,  is  determined  by  dividing  the  total  number  of  boundary 
pixels  (Nb)  by  the  desired  number  of  radius  vectors  (Nr)  i.e. 

Nh 

A  =  arclength  =  (4.3.5) 

r 

Thus  the  time  series  is  formed  from  the  radius  vector  lengths  correspond¬ 
ing  to  every  Ath  stored  boundary  pixel  where  the  length  from  pixel  to 
centroid  is  calculated  using  Equation  4.3.4.  This  method  is  simpler  and 
quicker  than  the  equi-angle  method  since  it  doesn't  need  to  check  each 
boundary  pixel  for  radius  vector  intersections. 


In  both  radius  vector  length  calculating  methods,  the  (x  ,y  ) 
position  of  the  boundary  intersection  pixel  of  the  maximum  radius  vector 
is  stored.  Also  the  length  of  the  maximum  radius  vector  is  stored. 

From  this  length  we  have  an  indication  of  the  size  of  shape.  Using  the 
position  we  can  calculate  the  orientation  of  the  maximum  radius  vector 


and  thereby  have  an  indication  of  the  orientation  of  entire  shape.  The 
orientation  of  maximum  radius  vector  in  the  Cartesian  coordinate  plane  is 


calculated  as 


orientation  (degrees)  =  Arctan 


(y  -y  1 
wmax  Jc ' 


—  xmax’xc 


The  negative  sign  is  due  to  the  origin  of  the  video  coordinate  system 
being  at  the  top  left  corner  of  the  scene. 

We  are  now  capable  of  forming  two  types  of  time  series  for  any  inner 
or  outer  object  boundary  encountered  in  the  scene.  Thus  it  is  possible 

to  calculate  the  AR  parameters  O,  ,9m2 . “mm’  a//e}  of  each  boundary 

at  any  model  order.  These  parameters  can  be  stored  in  vector  form  during 
the  acquisition  of  training  data  and  transformed  later,  or  they  can  be 
transformed  immediately  for  classification  purposes.  The  next  section 
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deals  with  the  extraction  of  additional  information  from  the  inner 
boundaries  of  an  object  and  how  that  information  can  be  used  in  the 
classification  process. 

4.4  Extraction  of  Additional  Information  from  the 
Inner  Boundaries  of  an  Object 

The  additional  information  provided  by  the  inner  boundaries  of  an 
object  can  be  used  to  hasten  the  classification  process  in  many  ways.  If 
the  number  of  inner  boundaries  or  holes  is  unique  i.e.  no  other  object 
has  that  number  of  holes,  then  the  classification  of  that  object  can  be 
based  solely  on  that  number.  Also  the  relative  position  of  the  edge 
used  in  conjunction  with  the  vector  of  transformed  AR  parameters  can 
provide  additional  information.  If  the  object  has  more  than  one  hole  then 
it  is  possible  to  calculate  the  AR  parameters  of  the  pr>gonal  shape  formed 
by  connecting  the  boundary  centroids  with  straight  lines.  In  this  work 
the  AR  parameters  of  the  polygon  were  not  calculated  but  the  number  of 
holes  and  the  relative  position  of  the  holes  inside  the  object  were  found. 
In  order  to  count  the  number  of  holes  belonging  to  an  object  it  is  first 
necessary  to  determine  which  holes  are  owned  by  which  object.  The  tech¬ 
nique  of  determining  hole  ownership  is  described  next. 

4.4.1.  Determination  of  Hole  Ownership 

Once  the  boundaries  have  been  detected  and  marked  as  in  Figure 

4.3.1,  hole  ownership  can  be  determined  using  the  following  technique. 

For  each  boundary  the  algorithm  is  as  follows: 

1)  Begin  at  the  start  pixel  and  proceed  to  read 
every  pixel  on  the  row  to  the  right  of  the 
start  pixel. 
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2)  Keep  track  of  any  occurrence  of  a  pixel  intensity 
which  is  neither  the  object  intensity  or  the  back¬ 
ground  intensity.  Once  a  pixel  intensity  denoting 
an  edge  is  detected,  ignore  any  occurrences  of  the 
same  intensity  if  they  happen  to  occur  immediately 
after  the  first  encounter. 

3)  Return  to  the  start  pixel  and  repeat  the  process  on 
the  same  row  but  now  going  towards  the  left. 

4)  Keep  a  separate  record  of  the  occurrences  of 
the  edge  pixels  as  in  Step  2). 

Once  both  right  and  left  sections  of  the  row  containing  the  start  pixel 
have  been  read,  it  is  then  possible  to  determine  which  boundary  the  edge 
in  question  is  inside.  If  the  same  boundary  intensity  is  detected  an  odd 
number  of  times  in  both  the  right  and  left  scans,  then  it  can  be  assumed 
that  the  edge  in  question  is  inside  the  boundary  corresponding  to  the 
intensity.  Figure  4.4.1  shows  an  example  of  how  this  method  can  be  used 
to  determine  that  edge  2  is  inside  boundary  1. 

1  crossing  detected  in  3  crossings  detected  in 

left  search  right  search 


Figure  4.4.1  Example  of  the  Method  for  Determining  Hole  Ownership 
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Once  the  inner  holes  of  an  object  have  been  identified  it  is  a  simple 
matter  to  count  them.  It  is  also  possible  to  order  them  according  to  the 
distance  between  their  centroids  and  the  centroid  of  the  outer  boundary. 

In  this  work  the  inner  edges  were  ordered  such  that  the  edges  whose 
centroids  were  furthest  away  from  the  outer  boundary  centroid,  were  put 
on  the  top  of  the  list.  When  the  inner  edges  have  been  detected,  counted, 
and  ordered,  it  is  then  possible  to  classify  an  object  simply  by  the  hole 
count  or  by  the  classification  of  each  one  of  its  edges. 

4.5  Description  of  System  Software 

The  algorithms  described  in  the  preceeding  pages  were  coded  into 
Fortran-77  for  the  VAX  11-730  computer.  The  two  main  programs  for  the 
pattern  recognition  system  are  called  TRAIN  and  CLASSIFY.  Program  TRAIN 
is  used  to  provide  the  classifier  with  the  training  data  for  each  ex¬ 
pected  object  edge.  Program  CLASSIFY  is  used  to  classify  objects  with 
one  or  more  edges  by  the  classification  of  each  edge  or  by  the  total 
number  of  inner  edges.  An  additional  program  called  MAKLIB  was  written 
to  create  a  library  file  (LIB. DAT)  of  AR  parameters  of  the  edges. 

This  library  file  increased  the  speed  of  system  tests  since  it  eliminated 
the  need  for  recalculating  the  AR  parameters  of  the  edges  for  each  test. 

A  complete  listing  of  TRAIN  and  CLASSIFY  and  the  subroutines  they  call  can 
be  found  in  the  Appendix  at  the  end  of  this  work. 

4.5.1  Program  TRAIN 

Program  TRAIN  is  an  interactive  program  which  trains  the  classifier 
with  the  data  of  the  edges  of  the  expected  objects.  This  program  pro¬ 
duces  a  TRAIN.DAT  file  which  contains  the  M+l  columns  of  the  trans- 


formation  matrix  and  also  the  average  of  the  transformed  samples  of  each 
class.  Multiple  edges  are  allowed  in  the  scenes,  however  all  edges  of  the 
objects  in  view  must  be  identical.  [See  Figure  4.5.1  for  the  general 
flowchart  of  program  TRAIN.]  The  program  begins  by  initializing  the  ITI 
hardware  and  then  prompts  the  user  for  the  scan  step,  the  AR  model 
order,  and  the  number  of  radius  vectors.  It  also  prompts  for  the  desired 
transformation  technique  which  can  be  either  the  basic  Karhunen-Loeve 
(KLSEL),  the  minimum  entropy  K-L  (MESEL),  or  the  optimal  prewhitened  K-L 
(PREKL)  coordinate  axes  generating  technique.  Once  the  scene  has  been 
viewed,  thresholded,  and  frozen,  it  is  scanned  in  order  to  find  the  start 
pixel  of  each  boundary.  As  soon  as  a  boundary  start  pixel  is  detected, 
the  entire  boundary  is  traced  and  the  centroid  is  calculated.  After  the 
entire  scene  has  been  scanned,  the  radius  vector  lengths  of  each  boundary 
are  determined  so  that  the  AR  parameters  can  be  calculated.  Once  the 
feature  vectors  are  formed  from  the  AR  parameters,  the  process  repeats 
and  another  scene  can  be  viewed.  When  enough  samples  of  each  class  have 
been  accumulated,  all  the  AR  parameter  feature  vectors  are  passed  to  the 
desired  subroutine  which  then  forms  the  average  covariance  matrix  and 
finds  the  particular  K-L  coordinate  axes.  After  the  K-L  axes  are  calculated 
and  the  transformation  matrix  is  formed,  the  original  feature  vectors  are 
transformed  and  the  average  transformed  vectors  are  found.  These  vectors 
are  stored  along  with  the  transformation  matrix  in  the  training  data 
file. 

4.5.2  Program  CLASSIFY 

Program  CLASSIFY  is  an  interactive  program  which  can  recognize  objects 


with  one  or  more  edges.  The  program  allows  many  different  objects  to  be 
present  in  the  scene.  It  also  calculates  the  rotation  of  the  objects.  The 
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Figure  4.5.1  General  Flowchart  of  Program  TRAIN 


program  starts  by  initializing  the  hardware  and  prompts  for  the  desired 
constants  as  in  Program  TRAIN.  It  also  prompts  for  the  desired  dimension 
of  the  average  transformed  feature  vectors.  Once  the  constants  have  been 
accepted,  the  program  scans  the  frame  buffer  and  turtles  around  all 
edges,  calculating  the  centroid  after  detecting  each  edge.  When  all  the 
edges  have  been  detected  and  also  marked,  the  program  then  determines  which 
holes  are  owned  by  which  objects.  The  inner  holes  of  the  objects  are 
counted  and  if  the  count  is  unique,  then  the  object  can  be  classified. 

The  edges  of  the  yet  unclassified  objects  now  need  to  be  classified  so 
the  radius  vector  lengths  are  determined  and  the  AR  parameters  are 
calculated.  Then  using  the  transformation  matrix  provided  by  the  training 
program,  the  AR  parameters  of  each  edge  are  transformed.  Now  the 
transformed  feature  vectors  can  be  classified  using  the  minimum-distance 
algorithm  as  discussed  in  Chapter  3.  Once  all  the  edges  of  the  object 
are  classified,  the  object  can  be  classified.  This  is  done  by  matching 
the  edge  classifications  and  the  number  of  holes  data  with  the  object 
data  contained  in  a  file  called  OBJECT.DAT.  This  data  file  is  created 
by  the  user  in  the  same  way  that  an  ordinary  text  file  is  created.  The 
object  data  is  typed  in  one  line  per  object.  For  example,  suppose  we 
have  an  object  with  label  3  which  has  2  holes  labelled  5  and  7.  If  the 
centroid  of  hole  7  is  further  away  from  the  outer  boundary  centroid  than 
hole  5,  then  we  would  type  in  the  sequence:  3275.  Once  all  the 
objects  have  been  classified,  the  program  is  ready  to  accept  another 
scene  full  of  objects  or  the  session  can  be  terminated.  See  Figure 
4.5.2  for  the  general  flowchart  of  Program  CLASSIFY. 
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Figure  4.5.2  General  Flowchart  of  Program  CLASSIFY 


CHAPTER  5 


NUMERICAL  EXAMPLE  AND  SYSTEM  TESTS  AND  RESULTS 

In  the  previous  chapters,  we  described  the  theory  and  implementation 
of  our  pattern  recognition  algorithms.  It  is  well  known  that  any 
theoretical  study  is  incomplete  without  examples  and  experimental  tests. 
Thus  this  chapter  is  divided  into  two  parts,  the  first  part  provides  a 
numerical  example  of  the  feature  selection,  transformation,  and 
classification  theory.  The  second  part  describes  the  series  of  tests 
and  their  results,  of  the  entire  recognition  system.  In  all  examples 
and  system  tests  the  image  threshold  was  set  at  120,  the  scan  step  was  8 
pixels,  and  the  desired  number  of  radius  vectors  was  64. 

5.1  Numerical  Example 

We  used  programs  TRAIN  and  CLASSIFY  to  provide  for  the  following 
examples  of  the  theories  and  algorithms  used  in  this  pattern  recognition 
system.  Two  examples  have  been  provided,  the  first  is  of  the  system  using 
the  equal  angle  boundary  sampling  method.  The  second  example  shows  the 
system  using  the  equal  arc  length  boundary  sampling  method.  The  shapes 
chosen  for  these  examples  are  the  familiar  S  (class  1)  and  N  (class  2) 
used  in  the  work  of  Dubois  [Dubois,  1984].  These  shapes  were  chosen  for 
two  reasons.  The  first  is  that  the  AR  parameters  of  both  classes  are 
very  similar  and  thus  they  show  the  need  for  a  type  of  transformation 

i 

technique.  Also  they  were  chosen  to  allow  quick  comparison  of  the  two 
different  pattern  recognition  systems.  The  actual  size  of  these  shapes 


are  considerably  larger  than  normal  text.  The  letter  S  is  approximately 
6  inches  tall  and  the  N  is  about  4-1/2  inches  tall. 

5.1.1  Example  of  the  System  Using  the  Equal  Angle  Boundary  Sampling  Method 

The  equal  angle  boundary  sampling  technique  described  in  Chapter  Two 
is  used  in  this  example.  Figure  5.1.1  shows  the  sampled  shapes  of  class  1 
and  class  2.  Note  that  portions  of  the  letters  are  not  sampled  as  would 
be  expected.  This  is  due  to  the  radius  vector  boundary  intersection 
algorithm  rejecting  radius  vector  measurements  if  they  occur  on  the  same 
slope  as  the  previous  radius  vector  length  measurement.  This  rejection 
was  done  to  prevent  the  occurrence  of  a  disproportionately  large  number  of 
similar  radius  vectors  measurements  if  the  slope  of  the  boundary  happened 
to  coincide  with  the  slope  of  a  radius  vector.  Below  the  shape  of  each 
letter  is  the  corresponding  time  series  formed  from  the  radius  vector 
lengths.  Notice  the  fast  vertical  changes.  These  unnatural  transitions 
are  the  direct  result  of  the  sampling  process  skipping  or  by-passing 
certain  sections  of  the  boundary.  These  quick  transitions  increase  the 
high  frequency  content  of  the  time  series  wave  shape  and  thus  it  can  be 
expected  that  the  resultant  AR  model  order  will  be  higher  for  this  type  of 
wave  shape  than  for  a  waveform  that  doesn't  have  these  sharp  transitions. 

Table  5.1.1  shows  the  original  training  set  data  comprising  10 
samples  of  each  class.  As  can  be  seen,  the  samples  varied  in  size  and 
rotation.  The  size  variation  can  be  observed  in  the  changes  of  the 
length  (in  pixels)  of  the  maximum  radius  vector  from  one  sample  to  the 
next.  Notice  also  that  the  changes  in  rotation  are  not  strict  multiples 
of  2-/N  .  This  was  done  to  mimic  a  realistic  setting  where  the  position 
and  rotation  of  an  object  is  generally  completely  arbitrary. 


Table 

5.1.1  Training 

Set  Data 

Equal  Angle  Sampling 

Method 

(AR  Model  Order: 

1) 

a//s" 

Length  of 

91 

Max  Rad  Vec. 

Rotation 

Class 

1 

0.97653 

1.73521 

183 

251.4 

0.97261 

1.8513 

187 

207.1 

0.96371 

1.82773 

188 

166.0 

0.97306 

1.80995 

187 

339.0 

0.96508 

1.79821 

193 

313.8 

0.96442 

1.86203 

196 

56.5 

0.95070 

2.10665 

203 

13.7 

0.96093 

1.81376 

204 

339.0 

0.97108 

1.67467 

214 

326.0 

0.91653 

1.92008 

209 

278.0 

Class 

2 

0.97338 

1.91497 

160 

236.0 

0.97753 

1.96800 

165 

193.8 

0.97996 

1.88319 

168 

159.2 

0.97825 

1.91980 

166 

326.0 

0.98253 

1.96457 

. 

174 

303.1 

0.97322 

1.93067 

177 

39.7 

0.97945 

1.94201 

182 

345.9 

0.97759 

1.92446 

190 

20.6 

0.97306 

1.93831 

186 

128.5 

0.97896 

1 .95099 

188 

213.5 

61 


The  AR  parameters  e-j  and  a//e  cluster  around  their  average  class 
value  as  can  be  seen  in  Figure  5.1.2a.  However  looking  at  Table  5.1.1, 
some  0-j  values  and  also  some  a//f  values  overlap  between  the  classes. 

Thus  it  can  be  seen  that  it  is  desirable  to  reduce  or  possibly  eliminate 
these  overlaps  by  using  a  transformation  technique.  Table  5.1.2  shows 
the  results  of  transforming  the  original  AR  parameters  using  the  basic 
K-L,  the  minimum  entropy  K-L  and  the  optimal  K-L  transformation  matrices 
as  described  in  Chapter  3.  The  plots  of  these  transformed  points  are 
shown  in  parts  b,c,d  of  Figure  5.1.2. 

To  understand  how  the  K-L  coordinate  axes  can  reduce  the  dimension 
of  the  feature  vector,  observe  how  the  transformation  rotates  the 
original  coordinate  axes  so  they  now  point  in  the  direction  of  maximum 
and  minimum  variance.  In  the  basic  K-L  transformed  system,  the  first 
axis  or  x  axis,  points  in  the  direction  of  maximum  within  class  variance. 
Looking  at  plot  (b)  in  Figure  5.1.2,  it  is  noticed  that  the  major  axis 
of  the  ellipsoid  formed  from  the  cluster  of  class  1  points  of  plot  (a) 
is  now  parallel  with  the  x  axis.  We  can  also  see  that  along  the  y  axis 
corresponding  to  the  direction  of  minimum  variance,  there  is  no  over¬ 
lapping  between  the  classes.  It  is  possible  to  draw  a  straight  line 
parallel  to  the  x  axis  between  the  class  clusters  such  that  the  classes 


do  not  overlap.  An  unlabelled  input  feature  vector  can  be  classified  by 
its  position  along  the  y  axis  alone.  Thus  the  dimension  of  the  feature 
space  has  been  reduced. 

Plot  (c)  in  Figure  5.1.2  shows  the  effect  of  the  minimum  entropy 
transformation.  Comparing  this  plot  with  plot  (b)  we  see  that  plot  (c) 
is  the  same  as  (b)  but  with  the  axes  reversed.  In  the  minimum  entropy 
transformed  system,  the  x  axis  points  in  the  direction  of  the  minimum 
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Table  5.1.2  Transformed  Training  Data 


Sample 

Basic 

K-L 

Min  Entropy  K-L 

Optimal 

K-L 

T1 

T2 

T1 

T2 

T1 

T2 

Class  1 

1 

1.685 

1.061 

1.061 

1.685 

263.35 

54.11 

2 

1.801 

1.063 

1 .063 

1.801 

269.21 

52.84 

3 

1.778 

1.053 

1.053 

1.778 

266.64 

52.41 

4 

1.760 

1.062 

1 .062 

1.760 

268.67 

53.23 

5 

1.748 

1 .053 

1.053 

1.748 

266.51 

52.76 

6 

1 .812 

1 .056 

1.056 

1.812 

267.35 

52.17 

7 

2.057 

1.054 

1.054 

2.057 

267.81 

49.09 

8 

1.764 

1.050 

1.050 

1.764 

265.73 

52.33 

9 

1.625 

1.053 

1 .053 

1.625 

266.06 

54.25 

10 

1 .870 

1 .056 

1 .056 

1 .870 

267.55 

51 .46 

Class  2 


1 

1.864 

1.067 

1.0671 

1.864 

270.40 

52.35 

2 

1.917 

1 .074 

1.074 

1.917 

272.25 

52.18 

3 

1.832 

1.072 

1.072 

1.832 

271.52 

53.08 

4 

1.869 

1.072 

1.072 

1.869 

271.68 

52.65 

5 

1.913 

1.079 

1.079 

1.913 

273.43 

52.56 

6 

1  .880 

1.067 

1.067 

1.880 

270.60 

52.20 

7 

1.891 

1.075 

1.075 

1.891 

272.32 

52.54 

8 

1.874 

1.072 

1.072 

1.874 

271.58 

52.56 

9 

1.888 

1 .068 

1 .068 

1 .888 

270.68 

52.12 

10 

1 .900 

1.075 

1.075 

1.900 

272.34 

52.43 

Figure  5.1.2  Plots  Showing  the  Effects  of  the  Different  Transformation 
Techniques:  (a)  Original  AR  Parameters,  (b)  Basic  K-L  Transformed 
Parameters,  (c)  Minimum  Entropy  K-L  Transformed  Parameters,  (d) 
Optimal  K-L  Transformed  Parameters 


within  class  variance.  This  is  to  be  expected  since  the  transformation 
matrix  is  formed  from  the  exact  same  eigenvectors  as  the  K-L  trans¬ 
formation  matrix,  but  reversed  in  order.  Since  this  reversal  of  eigen¬ 
vectors  does  not  enhance  the  basic  K-L  system,  the  minimum  entropy 
transformation  was  not  investigated  further.  However,  if  it  is  known  a 
priori  that  all  the  class  covariance  matrices  of  a  training  set  are  equal, 
then  it  is  shown  [Tou  and  Gonzalez,  1974]  that  the  minimum  entropy 
method  produces  better  results  than  the  basic  K-L  transformation. 

Plot  (d)  in  Figure  5.1.2.  shows  the  effect  of  the  optimal  K-L 
transformation.  The  x-axis  points  in  the  direction  of  maximum  between 
class  variance.  Now  the  class  means  fall  on  a  line  parallel  to  the  x- 
axis.  It  is  also  noted  that  the  intraset  distances  along  the  x-axis 
have  been  reduced  in  comparison  to  the  original  clusters.  Thus  any  un¬ 
labelled  input  vector  can  be  categorized  according  to  its  position  along 
the  x  axis. 

In  this  2  class  problem,  we  do  not  observe  a  distinct  advantage  of 
the  optimal  K-L  transformation.  However,  as  the  feature  space  increases 
and  includes  more  classes  the  ability  of  the  basic  K-L  transformation  to 
divide  the  feature  space  decreases.  This  is  because  the  K-L  trans¬ 
formation  is  based  on  the  within  class  variances  which  do  not  contain  a 
great  amount  of  discriminatory  information.  On  the  other  hand  the  optimal 
K-L  transformation  uses  the  prewhitened  between  class  variances  and  thus 
maximizes  the  use  of  the  discriminatory  information  contained  in  the 
class  means. 

Tables  5.1.3  and  5.1.4  show  the  data  of  the  classification  process. 
The  classifier  initially  retrieves  from  storage  the  transformation 
matrix,  which  was  used  to  transform  the  training  samples,  and  the 


Table  5.1.3  Classification  Data  Using  Basic  K-L  Transformation 
( AR  Model  Order:  1,  Transformed  Dimension:  2) 

INPUT  TRAINING  DATA 


Transformation  Matrix:  Mean  Transformed  Training  Samples: 


T1 

T2 

|ff.  04958 

-0.99877 

Class  1 

:  1.78981 

1.05600 

|0.99877 

0.04958 

Class  2 

:  1 .88286 

1.07206 

CLASSIFICATION  DATA 

a//s" 

Length  of 

01 

T1 

T2 

Max  Rad  Vec 

Rotation  Deci: 

Class  1 

-0.95892 

1 .85477 

1  .805 

1.050 

190 

66.8 

1 

-0.96278 

1.76381 

1.714 

1.050 

195 

27.0 

1 

-0.96952 

1.69998 

1 .650 

1.053 

195 

0.0 

1 

-0.97158 

1.88558 

1.835 

1.064 

192 

326.2 

1 

-0.97213 

1.75851 

1.708 

1.058 

209 

319.2 

1 

-0.97280 

1 .84680 

1.796 

1.063 

204 

122.9 

1 

-0.97630 

1.67335 

1.623 

1.058 

198 

246.8 

1 

-0.97096 

1 .80878 

1 .758 

1.059 

201 

213.3 

1 

-0.97648 

1.68312 

1.632 

1.059 

221 

213.7 

1 

-0.96265 

1.98290 

1.933 

1 .074 

219 

146.1 

2 

Class  2 

-0.97963 

1.92738 

1.876 

1.074 

180 

219.8 

2 

-0.98326 

1 .89282 

1 .842 

1 .076 

186 

172.7 

2 

-0.97366 

1.91872 

1 .868 

1 .068 

185 

7.0 

2 

-0.97945 

1.92982 

1.878 

1.074 

183 

314.0 

2 

-0.98020 

1 .94440 

1.894 

1.075 

198 

313.8 

2 

-0.97278 

1.941 15 

1 .891 

1.067 

197 

39.8 

2 

-0.98261 

1.95360 

1 .902 

1 .078 

200 

338.9 

2 

-0.97505 

1.93480 

1 .884 

1.070 

206 

13.9 

2 

-0.97879 

1 .91234 

1 .861 

1 .072 

208 

339.0 

2 

-0.97175 

1 . 91 964 

1.869 

1 .066 

205 

123.2 

2 

CONFUSION  MATRIX 
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191  95.0%  correct 


Table  5.1.4  Classification  Data  Using  Optimal  K-L  Transformation 


(AR  Model  Order:  1,  Transformed  Dimension:  2) 
INPUT  TRAINING  DATA 


Transformation  Matrix:  Mean  Transformed  Training  Samples: 


T1 

T2 

-246.90 

-70.75 

Class 

1  : 

267.388 

52.467 

_  15.70 

-  8.63 

Class 

2: 

271.682 

52.467 

CLASSIFICATION 

DATA 

Length  of 

el 

ot//8 

T1 

T2 

Max  Rad  Vec 

Rotation 

Dec i s i on 

Class  1 

-0.97576 

1 .79706 

269.13 

53.53 

183 

260.8 

1 

-0.97308 

1 .81335 

268.73 

53.20 

184 

231.3 

1 

-0.95263 

2.00491 

266.69 

50.10 

186 

186.9 

1 

-0.95786 

1.82341 

265.13 

52.04 

186 

139.8 

1 

-0.97406 

1.79275 

268.65 

53.45 

181 

94.4 

1 

-0.96344 

1.84091 

266.78 

52.28 

196 

71.7 

1 

-0.96077 

1.87224 

266.61 

51.82 

200 

39.8 

1 

-0.95672 

1.90406 

266.11 

51.26 

201 

345.8 

1 

-0.97174 

1.71986 

266.93 

53.91 

204 

326.0 

1 

-0.95385 

2.16893 

269.56 

48.77 

200 

288.2 

2 

Class  2 

-0.97945 

1.93141 

272.15 

52.63 

170 

134.2 

2 

-0.97635 

1.97708 

272.10 

52.02 

172 

200.7 

2 

-0.97666 

1 .88894 

270.80 

52.80 

173 

186.6 

2 

-0.97523 

1 .87643 

270.25 

52.81 

172 

152.5 

2 

-0.98230 

1.93227 

272.87 

52.83 

174 

0.0 

2 

-0.98010 

1.93268 

272.33 

52.67 

170 

297.8 

2 

-0.97795 

1.98704 

272.66 

52.05 

187 

33.7 

2 

-0.97538 

1.93191 

271.16 

52.34 

189 

346.0 

2 

-0.98182 

1.93898 

272.86 

52.73 

200 

288.2 

2 

-0.97969 

1.94883 

272.49 

52.50 

198 

123.2 

2 

CONFUSION  MATRIX 


1  2 

l 

1  j  9  1  95.0%  correct 

2  0  10 
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formed  training  sample  of  each  class.  The  table  shows  the  original  AR 
parameters  and  the  transformed  parameters  of  the  test  set.  The  maximum 
radius  vector  length  and  the  rotation  of  the  object  are  also  listed. 

Like  the  training  samples,  the  test  samples  are  random.  The  resultant 
classifier  decisions  are  shown  in  the  right  column.  For  both  sets  of 
classification  data,  the  basic  K-L  transformed  and  the  optimal  K-L 
transformed,  the  classifier  recognized  the  input  shape  correctly  for  all 
test  samples  except  one.  This  misclassification  occurred  in  both  sets 
for  a  sample  in  which  the  original  AR  parameters  vary  more  than  the 
rest.  It  is  suspected  that  this  difference  is  due  to  a  large  segment  of 
the  boundary  not  being  sampled.  It  is  expected  that  the  classification 
performance  at  this  low  model  order  (i.e.  M=l)  will  be  improved  by  the 
equal  arc  length  boundary  sampling  method.  This  will  be  investigated  next. 

5.1.2  Example  of  the  System  Using  the  Equal  Arc  Length  Sampling  Method 

The  equal  arc  length  boundary  sampling  method  as  described  in 
Chapter  Two  is  used  in  this  example.  Figure  5.1.3  shows  the  sampled 
shapes  of  class  1  and  class  2.  Below  each  sampled  shape  is  the  re¬ 
sultant  time  series.  As  can  be  observed  in  these  time  series  plots,  the 
sharp  vertical  transitions,  which  occurred  in  the  time  series  of  the 
previously  discussed  sampling  method,  do  not  occur.  The  waveforms  show 
the  gradual  changes  in  the  contour  of  the  boundary  as  it  is  sampled. 

Tables  5.1.5  and  5.1.6  show  the  training  and  classification  data.  For 
this  example,  only  the  optimal  K-L  transformation  was  used.  As  expected, 
the  AR  parameters  and  thus  the  transformed  parameters  cluster  closely 
about  their  respective  class  means.  As  a  result  we  see  a  perfect 


Table  5.1.5  Training  Set  DATA  Using  the  Equal  Arc  Length  Sampling 
Method  and  the  Optimal  K-L  Transformation 
(AR  Model :  1) 


Sample 

61 

a//e 

Length  of 

Max  Rad  Vec 

Rotation 

Transformed 

(Optimal 

T1 

1 

-0.97331 

2.26854 

Class  1 

202 

229.5 

-270.4 

2 

-0.97314 

2.31667 

203 

193.8 

-271.4 

3 

-0.97040 

2.34302 

200 

144.7 

-271.4 

4 

-0.97041 

2.46447 

204 

329.6 

-274.3 

5 

-0.97226 

2.38723 

219 

78.9 

-272.9 

6 

-0.97322 

2.29049 

217 

30.7 

-270.8 

7 

-0.97209 

2.35664 

216 

348.4 

-272.1 

8 

-0.96968 

2.44291 

220 

299.1 

-273.6 

9 

-0.97061 

2.43155 

225 

280.3 

-273.5 

10 

-0.97228 

2.30699 

183 

26.0 

-271 .0 

1 

-0.96430 

2.89333 

Class  2 

184 

111.5 

-282.9 

2 

-0.96269 

2.94088 

188 

180.0 

-283.6 

3 

-0.96187 

2.91875 

188 

137.5 

-282.9 

4 

-0.96357 

2.88612 

187 

335.1 

-282.6 

5 

-0.96345 

2.91985 

186 

78.7 

-283.3 

6 

-0.96316 

2.90567 

199 

270.6 

-282.9 

7 

-0.96178 

2.95423 

203 

336.4 

-283.7 

8 

-0.96132 

2.99844 

199 

296.0 

-284.7 

9 

-0.96314 

2.89153 

204 

276.2 

-282.6 

10 

-0.96405 

2.92276 

202 

113.8 

-283.5 

t 


* 


Samples 

K-L) 

T2 


1479.1 

1479.9 

1476.6 

1479.4 

1480.3 

1479.5 

1479.4 

1477.8 

1478.9 

1478.5 


1480.3 
1478.9 

1477.3 
1479.0 
1479.6 

1478.9 

1477.9 

1478.3 

1478.5 

1480.6 


Table  5.1.6  Classification  Data  Using  the  Equal  Arc  Length 
Sampling  Method  and  the  Optimal  K-L  Transformation 
(AR  Model  Order  1:  Transformed  Dimension:  1) 


INPUT  TRAINING  DATA 


Transformation  Matrix:  Mean  Transformed  Training  Samples: 


T1 

T2 

223.52 

-1466.18 

Class  1 : 

-272.15 

1478.9 

^23.27 

22.96 

Class  2: 

-283.27 

1478.9 

CLASSIFICATION 

1  DATA 

Length  of 

A 

V1 

a//g 

T1 

Max  Rad  Vec 

Rotation 

Decision 

Class  1 

-0.97266 

2.35135 

-272.14 

186 

260.9 

1 

-0.97328 

2.28844 

-270.81 

188 

231.0 

1 

-0.97302 

2.29317 

-270.86 

190 

192.4 

1 

-0.97031 

2.40888 

-272.95 

186 

146.9 

1 

-0.97180 

2.33380 

-271.54 

197 

147.7 

1 

-0.97035 

2.44807 

-273.87 

196 

317.8 

1 

-0.97194 

2.36741 

-272.35 

198 

77.9 

1 

-0.97324 

2.29135 

-270.87 

200 

26.4 

1 

-0.97202 

2.28805 

-270.52 

208 

27.5 

1 

-0.97262 

2.30752 

-271.11 

206 

346.5 

1 

-0.96444 

2.88873 

Class  2 

-282.81 

170 

289.9 

2 

-0.96061 

2.93651 

-283.06 

177 

8.9 

2 

-0.96100 

2.93754 

-283.17 

176 

321.7 

2 

-0.96298 

2.89563 

-282.64 

170 

276.5 

2 

-0.96398 

2.91838 

-283.39 

169 

242.3 

2 

-0.96461 

2.88745 

-282.82 

179 

114.0 

2 

-0.96244 

2.91625 

-283.00 

183 

199.2 

2 

-0.96240 

2 . 94407 

-283.64 

184 

142.2 

2 

-0.96264 

2.91241 

-282.96 

191 

12.8 

2 

-0.96283 

2.92158 

-283.21 

190 

295.6 

2 

CONFUSION  MATRIX 


1  2 


1 

2 


10 

0 


0 

10 


100%  correct 


classification  performance  using  only  one  dimension  of  the  transformed 
space.  This  example  shows  how  the  equal  arc  length  sampling  method 
produces  excellent  results  at  a  low  AR  model  order. 


5.2  System  Tests  and  Results 

We  used  programs  TRAIN,  CLASSIFY  and  MAKLIB  to  test  the  performance 
of  the  pattern  recognition  system.  We  tested  the  system  with  different 
shape  categories  and  with  varying  numbers  of  shapes  in  each  category.  In 
all  the  tests  the  training  data  for  all  the  desired  model  orders  was 
acquired  from  a  library  data  file  containing  the  AR  parameters  for  model 
orders  1-10  for  the  samples  of  the  particular  shape  category.  Thus  for 
each  shape  category  the  training  data  for  model  orders  1  through  10  was 
formed  from  the  same  set  of  samples. 

For  each  test  we  have  provided  a  table  of  results  showing  the  overall 
percent  correct  classifications  of  the  shapes  category  at  a  certain 
transformed  dimension  and  particular  model  order.  The  percent  calculation 
value  was  calculated  as  the  average  of  the  class  percent  correct  classi¬ 
fications.  This  method  of  calculating  percentages  downgrades  the  classifi 
more  if  the  misclassifications  occur  for  one  particular  class  than  if  the 
misclassifications  occurred  randomly  among  all  the  classes. 

All  the  shapes  used  in  the  tests  were  cut  out  of  white  paper  and 
viewed  against  a  black  background.  This  was  done  to  eliminate  the 
effects  of  shadows,  reflections,  and  texture  on  the  performance  of  this 
pattern  recognition  system.  In  all  the  tests  the  objects  were  translated 
and  rotated  randomly  in  the  scene.  Also  the  camera  was  lowered  and 
elevated  randomly  between  samples  to  change  size  of  the  viewed  object. 

In  all  the  tests,  except  the  multi -edge  object  and  timed  tests,  the 


classification  data  was  provided  by  the  same-  library  file  as  was  used  in 
the  training  process.  This  was  done  to  speed  up  the  tests  since  the  AR 
parameters  could  simply  be  pulled  out  of  storage  instead  of  calculated. 

We  felt  that  the  system  would  still  give  accurate  results  since  the  only 
data  the  classifier  uses  in  its  decision  making  process  is  the  mean  class 
transformed  training  samples  where  the  means  are  calculated  over  25 
samples. 

5.2.1  Industrial  Shapes  Test 

This  was  the  first  test  done  and  it  was  used  to  compare  the  basic 
K-L  vs  the  optimal  K-L  transformation  techniques.  We  chose  this  set  to 
again  allow  comparison  of  our  system  with  the  system  of  Dubois.  The  8 
industrial  shapes  (see  Figure  5.2.1)  in  this  test  are  scaled  down  versions 
of  the  industrial  shapes  tested  by  Dubois.  These  shapes  were  originally 
collected  from  the  literature  concerning  industrial  shape  pattern 
recognition. 

Table  5.2.1  shows  the  classification  results  of  the  system  using 
both  the  basic  K-L  transformation  and  the  optimal  K-L  transformation.  As 
can  be  observed  the  basic  K-L  transformed  system  produces  classification 
rates  of  91.5  to  98.5  percent  at  the  maximum  transformed  dimension  for 
model  orders  1  through  10.  A  significant  improvement  is  seen  in  the 
classification  results  of  the  optimal  K-L  transformed  system.  At  model 
orders  3  through  10  we  see  perfect  classification  rates  and  for  model 
orders  5,7,  and  9  these  perfect  rates  occur  when  using  only  two  dimensions 
of  the  transformed  feature  vector  space.  Generally  the  misclassifications 
occurred  when  the  classifier  confused  objects  1  and  8  and  objects  4  and  5. 


Table  5.2.1  Results  of  the  Industrial  Shapes  Test  Using  the  Basic  K-L 
and  the  Optimal  K-L  Transformation  Matrices 


(25  Test  Samples  =  25  Training  Samples) 
%  Correct  Classification 


Transformed  Dimension 


Model 

Order 

1 

2 

3 

4 

5 

6 

7 

8 

9 

10 

BASIC 

K-L 

1 

1 

j 

1 

96.0 

96.0 

2 

91.5 

92.0 

97.0 

3 

90.5 

91.0 

93.5 

93.5 

4 

94.5 

91 .0 

91 .0 

91.5 

91.5 

5 

89.0 

90.0 

95.5 

96.0 

96.0 

96.0 

6 

49.0 

87.5 

88.5 

91.5 

91.5 

91.5 

91.5 

7 

62.5 

89.5 

89.5 

96.5 

96.5 

98.0 

98.5 

98.5 

I 

8 

58.5 

78.0 

80.5 

89.0 

93.5 

97.5 

97.5 

97.5 

98.0 

9 

83.0 

84.5 

86.0 

84.0 

92.0 

94.0 

96.0 

96.5 

97.0 

97.0 

OPTIMAL 

K-L 

1 

85.0 

94.5 

2 

87.5 

98.0 

99.0 

3 

86.5 

95.5 

100.0 

100.0 

4 

80.0 

94.5 

100.0 

100.0 

5 

77.0 

100.0 

100.0 

100.0 

6 

84.5 

99.0 

100.0 

100.0 

7 

96.0 

100.0 

100.0 

100.0 

8 

88.0 

99.0 

98.5 

99.0 

9 

91 .0 

100.0 

100.0 

100.0 

10 

90.5 

98.0 

99.5 

100.0 

100.0 

100.0 

100.0 

100.0 

100.0 

100.0 

100.0 

100.0 

100.0 

100.0 

100.0 

100.0 

100.0 

100.0 

100.0 

100.0 

100.0 

100.0 

100.0 

100.0 

100.0 

100.0 

100.0 

100.0 

100.0 

100.0 

100.0 

This  test  checked  the  system  performance  on  a  category  of  4  military 
shapes  as  shown  in  Figure  5.2.2.  One  aspect  of  this  test  was  to  demonstrate 
how  this  system  could  be  used  for  military  target  classification.  Also, 
since  all  the  silhouettes  are  long  and  pointed,  another  aspect  was  to 
see  how  the  system  would  perform  in  the  classification  of  fairly  similar 
shapes.  The  optimal  K-L  transformation  matrix  was  used  and  both  boundary 
sampling  methods  were  implemented  in  this  test. 

The  test  results  in  Table  5.2.2  show  that  the  classifier  performed 
perfectly  using  the  equal  angle  sampling  method  and  1  dimension  of  the 
transformed  feature  space.  The  classifier  performed  slightly  less  well 
for  the  equal  arc  length  sampling  method.  The  misclassifications 
occurred  when  the  classifier  confused  shapes  1  and  3.  Note  that  shape  1  is 
the  smallest  of  the  shapes  in  this  category. 

5.2.3  Geometric  Shapes  Test 

Many  two  dimensional  representations  of  objects  are  fairly  geo¬ 
metrical  in  shape,  i.e.  the  shapes  are  symmetric  and  composed  of  lines 
and  curves.  Thus  a  test  was  performed  on  a  set  of  8  geometric  shapes 
(see  Figure  5.2.3).  This  test  also  provided  the  individual  edge 
classification  results  which  could  then  be  compared  with  the  results  of 
the  next  test  where  the  edges  are  combined  to  form  multiple  edge  shapes. 

As  in  the  military  shapes  test,  the  optimal  K-L  transformation  matrix 
and  both  boundary  sampling  methods  were  implemented. 

Table  5.2.3  shows  the  classification  results.  Note  that  the 


classification  rate  ranges  from  93.5%  to  99.5%  for  the  maximum  transformed 
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Table  5.2.2  Results  of  the  Military  Shapes  Test  Using  the  Optimal 
K-L  Transformation  and  Both  Boundary  Sampling  Methods 
(25  Test  Samples  =  25  Training  Samples) 

%  Correct  Classification 


Model  Transformed  Dimension 

Order  1  2 

EQUAL  ANGLE  SAMPLING  METHOD 

1  100.0  100.0 

2  100.0  100.0 


EQUAL  ARC  LENGTH  SAMPLING  METHOD 
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Table  5.2.3  Results  of  the  Geometric  Shapes  Test  Using  the  Optimal  K-L 
Transformation  and  Both  Boundary  Sampling  Methods 


(25  Test  Samples  =  25  Training  Samples) 
%  Correct  Classification 


Transformed  Dimension 


Model 

Order 

1 

2 

3 

4 

5 

6 

7 

8 

1 

95.5 

98.0 

2 

53.5 

94.5 

98.5 

3 

79.5 

96.0 

97.0 

93.5 

4 

51.5 

92.0 

98.5 

98.5 

95.0 

5 

68.0 

99.0 

99.5 

98.0 

98.5 

98.5 

6 

58.5 

93.5 

98.5 

99.0 

99.0 

99.0 

99.5 

7 

77.5 

95.0 

98.0 

99.0 

98.5 

99.0 

99.0 

99.0 

8 

80.5 

95.0 

97.0 

98.5 

98.5 

99.0 

99.0 

99.0 

EQUAL  ARC  LENGTH 

SAMPLING  METHOD 

1 

89.0 

98.5 

2 

72.5 

91.5 

93.0 

3 

58.5 

82.5 

83.5 

84.5 

4 

71.0 

86.5 

87.5 

88.5 

89.0 

5 

83.0 

89.5 

85.5 

86.5 

87.0  87.0 

9 


99.0 


dimensions  of  the  system  using  the  equal  angle  sampling  method.  Note 
also  that  the  best  recognition  rate  for  the  equal  arc  length  sampling 
method  occurs  at  the  lowest  model  order.  It  is  evident  upon  looking  at 
the  top  half  of  Table  5.2.3  that  2  samples  (from  class  5)  of  the  library 
data  file  had  AR  parameters  that  deviated  too  far  from  the  class  mean  and 
thus  were  consistently  misclassified.  If  the  test  had  used  all  random, 
new  samples,  we  would  probably  see  an  occasional  perfect  classification 
score.  However,  the  table  does  show  how  the  AR  parameters  can  deviate 
so  far  from  the  class  mean  that  even  the  best  transformation  technique 
cannot  completely  eliminate  the  effect  of  this  variance. 

5.2.4  Multi -Edge  Shape  Test 

This  multi-edge  shape  test  is  simply  an  extension  of  the  geometric 
shapes  test.  It  demonstrates  how  the  pattern  recognition  system  can  be 
used  to  recognize  shapes  with  multiple  edges.  For  this  test  we  produced 
12  shapes  with  multiple  edges  by  cutting  geometrical  shapes  out  of  the 
geometric  shapes  of  the  previous  section.  See  Figure  5.2.4  for  the  re¬ 
presentative  shapes.  The  training  data  for  this  test  was  formed  from  the 
same  samples  that  were  used  in  the  previous  test.  During  classification 
however,  the  test  data  was  formed  from  random,  new  samples.  Thus  this  is 
the  first  test  where  the  test  and  training  data  samples  are  mutually 
exclusive. 

The  confusion  matrix  of  one  run  of  this  test  at  model  order  5  and 
transformed  dimension  2  is  shown  in  Figure  5.2.5.  The  confusion  matrix 
shows  the  number  of  each  of  the  classifier  decisions  for  a  shape,  where 
each  shape  was  submitteJ  to  the  classifier  20  times.  It  is  observed  that 
some  samples  were  not  classified  at  all,  especially  samples  of  shape  1. 


82 

Figure  5.2.5  Confusion  Matrix  of  the  Multi-Edge  Shape  Test 
(AR  Model  Order:  5,  Transformed  Dimension:  2) 

(20  Test  Samples;  25  Training  Samples) 
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This  problem  was  caused  by  the  individual  edge  classifications  not  matching 
the  actual  object  vector  edge  classifications.  This  apparently  high 
rate  of  misclassification  is  due  mostly  to  the  greatly  reduced  size  of 
the  interior  edges.  When  the  size  of  an  edge  becomes  very  small  the  pixel 
quantization  problem  becomes  increasingly  significant.  The  resultant 
increased  noise  decreases  the  overall  signal  to  noise  ratio  which  ulti¬ 
mately  degrades  system  performance.  The  most  common  effect  of  the 
lowered  signal  to  noise  ratio  was  the  classifier  confusion  of  edge  6  with 
5  in  object  1. 

5.2.5  Large  Number  of  Shapes  Test 

In  this  test  we  investigated  the  performance  of  the  system  when  the 
classifier  was  given  a  large  training  set.  For  this  test  the  sample  data 
files  of  the  industrial,  military  and  geometric  shapes  were  combined  to 
form  a  single  20  shape  data  file.  The  classifier  was  trained  and  tested 
with  the  samples  of  this  data  file  using  the  optimal  K-L  transformation 
technique.  All  the  AR  parameters  in  the  data  file  were  formed  from  the 
time  series  produced  by  the  equal  angle  boundary  sampling  method. 

The  results  of  this  test  are  shown  in  Table  5.2.4.  The  individual 
runs  of  this  test  started  at  model  order  5  since  it  was  desired  that  a 
recognition  rate  of  at  least  99%  occur  for  the  individual  categories  be¬ 
fore  the  combination  was  tested.  As  can  be  seen,  the  classification  per¬ 
formance  degrades  somewhat  in  comparison  to  the  individual  tests.  This 
was  expected  since  the  feature  vector  space  becomes  more  crowded  as  more 
shapes  are  included  in  the  training  set. 


Correct  Classification 


Model 

Order 


Transformed  Dimension 


1 

2 

3 

4 

5 

6 

7 

53.2 

88.0 

93.8 

95.2 

95.4 

95.6 

57.8 

89.2 

91.4 

95.0 

95.4 

96.2 

97.0 

56.2 

90.6 

94.0 

95.4 

95.8 

97.0 

96.6 

44.2 

83.8 

96.2 

90.8 

92.6 

94.2 

94.2 
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5.2.6  Classification  Speed  • 

In  the  previous  sections  we  have  investigated  the  accuracy  of  the 
pattern  recognition  system.  Another  equally  important  aspect  of  a 
pattern  recognition  system  is  classification  speed.  Usually  the  classifier 
is  trained  off-line,  since  it  is  an  interactive  process,  so  training  speed 
is  desirable  but  not  critical.  What  is  critical  and  what  also  determines 
the  possible  applications  of  a  pattern  recognition  system  is  the  speed  of 
classification  combined  with  the  accuracy. 

In  the  process  of  classifying  a  single  edge,  the  edge  is  first 
detected  and  traced.  Then  the  radius  vectors  are  calculated  and  the  AR 
parameters  are  recursively  estimated.  The  AR  parameters  are  then  trans¬ 
formed  and  categorized  according  to  the  nearest  class  mean  in  the  feature 
space.  For  multiple  edge  objects  this  process  is  repeated  for  each  edge. 
Then  all  the  individual  edge  classifications  are  ordered  and  matched  to 
the  correct  object  vector  as  described  in  Chapter  Four.  If  the  multi  - 
edge  object  has  a  unique  number  of  inner  boundaries,  then  it  is  not 
necessary  to  calculate  the  AR  parameters  thus  reducing  the  classification 
time. 

The  most  time  consuming  tasks  of  the  classification  process,  using  the 
equal  angle  boundary  sampling  method,  are  ordered  below  with  the  most  time 
consuming  task  first: 

1)  scanning  the  entire  frame  and  tracing  the 
boundaries 

2)  radius  vector  length  calculation 

3)  AR  parameter  estimation 


The  rest  of  the  processes:  determining  hole  ownership,  ordering  the 
boundaries,  and  actual  classification  do  not  contribute  a  significant 
amount  of  time  to  the  total  classification  process.  Each  of  the  above 
tasks  are  accomplished  within  a  range  of  times  depending  on  the  factors 
as  listed  next  to  each  task  number  below: 

1)  scan  step  size,  number  of  edges,  total  length 
of  each  edge 

2)  number  of  desired  radius  vectors,  length  of  the 
edge  (equal  angle  sampling  method) 

3)  number  of  resultant  radius  vectors,  and  AR 
model  order. 

The  observed  times  for  scanning  the  scene  and  turtling  the  boundaries 
range  from  0.25  seconds,  for  no  objects  in  the  scene  (scan  step  =8),  to 
greater  than  2.2  seconds  for  a  scene  with  many  objects  with  multiple  edges. 
Table  2.25  shows  the  radius  vector  and  AR  parameter  calculation  times  for 
the  class  1  (S)  shape  of  the  numerical  example  in  the  first  section  of  this 
chapter.  The  length  of  the  boundary  for  this  shape  was  1308  pixels  and  the 
scan  and  turtle  time  was  approximately  0.83  seconds. 

Average  total  times  for  the  classifications  of  the  large  shapes  of 
the  numerical  example  are  1.32  seconds  using  the  equal  angle  sampling 
method  and  1.01  seconds  for  the  equal  arc  length  sampling  method.  For 
smaller  shapes  the  total  classification  time  for  a  single  edge  is  less 
than  one  second. 

In  this  chapter  we  used  the  programs  described  in  Chapter  Four  to 
demonstrate  the  system  theories  and  the  system  accuracy.  In  the  first 
section  we  present  two  examples  of  the  feature  selection  and  feature 
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Table  5.2.5  Processor  Times  for  Calculating  the  Radius 
Vectors  and  the  AR  Parameters 
(Fractions  of  a  Second) 
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Model 

Order 


Equal  Angle 
Sampling  Method 

.3398 

.3398 

.3515 

.3789 

.3789 

.3984 

.4023 

.4102 

.4180 

.4336 


Equal  Arc  Length 
Sampling  Method 

.0391 

.0508 

.0560 

.0703 

.0781 

.0781 

.0898 

.1094 

.1211 

.1211 
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vector  transformation  methods.  The  second  section  describes  a  series  of 
system  tests  and  the  resultant  classifier  performance.  In  the  last  part 
of  the  tests  and  results  section  we  summarize  the  time  requirements  of 
this  pattern  recognition  system. 


CHAPTER  6 


CONCLUSIONS  AND  SUGGESTED  IMPROVEMENTS 

In  this  thesis  we  have  developed  a  pattern  recognition  system 
capable  of  classifying  two  dimensional  shapes  with  many  edges.  Also  the 
system  is  capable  of  classifying  many  shapes  in  a  single  scene.  Since  the 
problem  presented  by  the  existence  of  many  edges  in  a  shape  can  be  con¬ 
sidered  as  a  simple  extension  of  the  single  edge  problem,  we  proceeded  to 
concentrate  on  improving  the  single  edge  classification  process.  In  so 
doing,  we  used  the  work  of  Dubois  [Dubois,  1984]  as  a  starting  point  and 
continued  by  improving  some  of  the  existing  algorithms  and  developing  a 
more  efficient  classification  scheme.  By  applying  the  theory  of  the 
Karhunen-Loeve  transformation  onto  the  feature  vectors  formed  from  the  AR 
model  parameters  of  an  objects  boundary,  we  were  able  to  reduce  the 
dimension  of  the  feature  vector  space  and  thus  further  reduce  training 
data  storage  space  and  also  classification  time. 

We  have  split  the  contents  of  this  chapter  into  two  sections.  The 
first  section  provides  comments  on  the  work  just  completed.  The  second 
section  suggests  future  work  that  can  be  done  to  enhance  the  performance 
of  our  system. 


6. 1  Concluding  Remarks 

On  the  feature  selection  techniques: 

The  recursive  estimation  of  the  AR  parameters  speeded  up 
classification  time  without  degrading  the  system  performance. 
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It  is  noted  that  most  of  the  time  spent  in  the  production  of 
the  AR  parameters  of  a  given  time  series  is  not  in  the  recursion 
itself,  but  in  calculating  the  correlation  coefficients  used  by 
the  recursive  equations. 

The  recursive  algorithm  did  not  alter  the  invariance  properties 
of  the  AR  parameters. 

The  performance  of  the  AR  model,  and  thus  of  the  classifier,  is 
dependent  on  the  size  of  the  shape.  When  the  size  of  the  shape 
as  viewed  on  the  monitor  decreases,  the  pixel  quantization  error 
increases.  Thus  the  noise  on  the  samples  of  the  time  series 
increases  and  the  AR  parameters  will  deviate  further  from  the 
class  means. 

The  equal  angle  radius  vector  length  calculation  method  of 
Dubois  was  speeded  up.  The  average  time  was  previously  1.68 
seconds,  this  can  be  compared  with  the  average  time  of  .33 
seconds  on  our  system.  However,  this  decreased  time  is  partly 
at  the  expense  of  an  increased  temporary  storage  requirement 
for  the  boundary  pixels. 

The  equal  arc  length  boundary  sampling  technique  can  be 
successfully  implemented  at  very  low  AR  model  orders.  This 
technique  results  in  an  average  .31  second  reduction  in 
classification  time  compared  with  the  time  of  the  equal  angle 
sampling  method. 
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The  equal  angle  boundary  sampling  method  consistently  provided 
better  results  than  the  equal  arc  length  sampling  method  at  AR 
model  orders  greater  than  1 . 

On  the  transformation  techniques: 

The  basic  and  the  optimal  Karhunen-Loeve  transformation 
techniques  are  effective  in  decorrelating  the  components  of 
the  feature  vector.  By  removing  the  correlated  components 
corresponding  to  the  smallest  eigenvalues  of  the  K-L  covariance 
matrices,  the  original  feature  vectors  are  approximated  by 
vectors  of  lower  dimension. 

The  optimal  K-L  transformation  technique  proves  to  be  better 
than  the  basic  K-L  transformation  when  there  are  more  than  2 
classes  in  the  feature  space.  This  is  because  the  optimal 
transformation  is  based  on  a  between  class  scatter  matrix 
whereas  the  basic  K-L  transformation  is  based  on  a  within 
class  scatter  matrix.  In  the  optimal  K-L  system  the  scatter 
matrix  is  prewhitened  to  remove  the  individual  sample  vari¬ 
ances.  This  is  done  to  provide  maximum  use  of  the  dis¬ 
criminatory  information  contained  in  the  class  mean  vectors. 

In  the  industrial  shapes  test  the  optimal  K-L  transformation 
technique  proved  to  be  more  effective  in  separating  the  feature 
space  than  the  rotated  coordinate  system  of  Dubois.  For  this 
test  the  shapes,  which  were  scaled  down  versions  of  the  shapes 
used  by  Dubois,  were  correctly  classified  at  AR  model  orders  3 
through  10.  Also  model  orders  5,7,9  required  only  2  dimensions 


of  the  transformed  space  to  achieve  correct  classification. 

This  can  be  compared  with  the  results  of  a  similar  test  by 
Dubois  where  correct  classification  was  achieved  only  at  model 
order  5  and  maximum  dimension  of  the  feature  space. 

As  in  the  rotated  coordinate  system  in  the  work  of  Dubois,  the 
K-L  transformation  matrices  are  based  on  the  eigenvectors  of  the 
covariance  matrices  of  the  feature  vectors.  These  covariance 
matrices  become  singular  whenever  the  number  of  training 
samples  is  less  than  M+2,  where  M  is  the  model  order.  Since  we 
consistently  used  25  training  samples  for  all  our  system  tests 
we  did  not  encounter  any  singular  matrices.  However  we  would 
like  to  note  that  this  singular  matrix  problem  can  very  well 
occur,  especially  at  high  model  orders. 

We  did  not  concern  ourselves  with  finding  the  optimal  model 
order  for  the  tests  since  the  shapes  used  in  the  tests  were  so 
varied.  Generally  the  best  results  were  obtained  with  model 
orders  3  through  10  and  with  model  order  1  for  the  equal  arc 
length  boundary  sampling  method. 

The  permanent  storage  requirements  for  the  trained  classifier 
in  this  system  consists  of  an  (m+1)  x  (m+1)  dimensioned  trans¬ 
formation  matrix  and  c  (m+1)  transformed  class  mean  feature 
vectors.  As  we  saw  in  the  examples  of  Sections  5.1.1  and  5.1.2 
where  c=2  and  m=l ,  the  storage  requirement  consisted  of  8 
values.  This  aspect  of  th's  pattern  recognition  system  makes 
it  ideal  for  the  classification  of  a  large  set  of  shapes  since 
the  storage  requirements  for  each  shape  is  minimal. 
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On  the  original  work  in  this  thesis: 

A  boundary  marking  and  bookkeeping  system  was  devised  to  allow 
multiple  shapes  with  multiple  edges  in  the  scene. 

A  simple  method  of  determining  hole  ownership  was  devised. 

The  theory  of  the  Karhunen-Loeve  expansion  was  applied  to  the 
vectors  formed  from  the  parameters  of  the  AR  model  of  the  edges 
of  a  shape. 

A  special  purpose  program  (MAKLIB)  was  written  to  collect  the 
AR  parameters  at  model  orders  1  through  10  for  any  desired 
edge.  This  collection  of  AR  parameters  was  used  for  classifier 
training  and  testing. 

Finally,  a  complete  set  of  software  (TRAIN  and  CLASSIFY)  was 
developed  into  a  pattern  recognition  system  which  not  only 
classifies  shapes  but  also  produces  size  and  rotation  indicators 
based  on  the  length  and  position  of  the  maximum  radius  vector. 

As  initially  stated,  we  set  out  to  classify  objects  with  multiple  edges. 

We  conclude  that  we  were  successful. 

6.2  Suggested  Improvements 

The  two  most  important  aspects  of  a  pattern  recognition  system  are 
speed  and  accuracy.  Hence  it  is  the  goal  of  this  section  to  present 
suggestions  which  will  increase  the  classification  speed  and  accuracy. 

To  increase  the  classification  speed: 

Implement  the  turtle  algorithm  in  assembly  language. 
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For  multiple  edge  objects  -  find  the  AR  parameters  of  the 
shape  formed  by  connecting  the  centroids  with  straight  line 
segments.  This  will  avoid  the  time  needed  to  classify  each 
individual  edge. 

Segment  the  scene  before  scanning  to  quickly  determine  the 
approximate  location  of  an  object. 

To  increase  classification  accuracy: 

Improve  the  equal  angle  boundary  sampling  technique  so  that  it 
does  not  skip  over  sections  of  the  boundary. 

Determine  the  optimum  AR  model  order  by  using  Equation  2.2. 15e 
as  described  in  Chapter  Two  and  then  use  the  value  as 
additional  information  in  the  classification  process. 

Also  to  increase  the  versatility  of  the  system: 

Modify  program  TRAIN  so  that  it  can  train  on  entire  objects 
without  having  to  isolate  the  edges. 

Investigate  different  thresholding  techniques  to  deal  with 
the  reflection,  shadow,  and  texture  problems  of  real  three 
dimensional  objects. 
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CAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA 

C 

C  FILE  TRAIN. FOR  LAST  REVISION  JULY  1,1985 

C  AUTHOR:  RUTH  D.  KENNETT 


THIS  FILE  CONTAINS: 


PROGRAM  TRAIN 


TO  LINK  TYPE: 


LINK  TRAIN,COLLECT,UORKl,RADl,PARAM,REGIS,KLSEL,MESEL,PREKL,EIGV,- 
USPKD,SYS*LIBRARY:VIDEOLIB 


caaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaa 
C  PROGRAM  TRAIN 

caaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaa 

c 

C  THIS  PROGRAM  TRAINS  THE  CLASSIFIER  BY  INTERACTIVELY  TAKING  SAMPLES 

C  OF  THE  CLASSES  REPRESENTING  THE  EDGES  OF  THE  EXPECTED  OBJECTS. 

C  ONCE  THE  USER  HAS  THE  DESIRED  NUMBER  OF  SAMPLES  OF  EACH  CLASS,  THE 

C  SAMPLES  ARE  THEN  SENT  TO  ONE  OF  THE  DESIRED  TRANSFORMATION  ROUT- 

C  INES:  KLSEL.MESEL.PREKL. 

C 


INTEGER  DEV, IN, OUT, STEP, INRAD, M 

INTEGER  D  IM, NSAMP, NCLASS, CLASS, BNDCT, SAMP 

INTEGER  PCLASS , TOTCLS 

INTEGER  INCODE, IRCODE 

INTEGER  NEULUN,OLDLUN 

INTEGER  XPIX(2500,10),YPIX(2500,10> 

INTEGER  BOUND (10,10) 


FAD< -2: 250, 10), ORIENT  (10) 

X  ( 1 1 , 25 ,20 ) , THETA ( 11,10) 

TRANS ( 11, 11), 1(25,11, 20) 

TA( 1 1 , 20) 

LIB(20,2S,10,11),ALL(10,11,10) 


CHARACTER*!  INST 


COMMON  DEV,  IN, OUT, STEP, INRAD, M 
DATA  IN, OUT/5, 5/ 

DATA  TRANS/121A0.0/ 


NEULUN  -  1 
QLDLUN  =  2 

OPEN(UNII=NEULUN,FIL£=  'TRAIN.DAT',  STATUS= ' NEU ' ) 
OPEN(UNIT=OLDLUN,FILE='LIB.DAT' , STATUS" ' OLD ' > 


SETUP  THE  IMAGING  HARDWARE  AND  PROMPT  FOR  DESIRED  PARAMETERS 


CALL  SETUP 
CALL  THRESH 


URITE(OUT, A)  'TYPE  1  -  TO  TRAIN  WITH  NEU  SAMPLES' 

UR  ITE ( OUT , A  >  'TYPE  2  -  TO  TRAIN  WITH  SAMPLES  FROM  LIB. DAT' 


READ( IN , 2 )  INCODE 
FORMAT ( 12) 


UR  ITE ( OUT , A )  'TYPE  1  -  FOR  MESEL  (MINIMUM  ENTROPY)' 


mm 


URITE(OUI,A>  'TYPE  2  -  EQR  KLSEL  ( KARHUNEN-LOEVE > ' 

UR ITE ( OUT , A )  'TYPE  3  -  FOR  PREKL  (OPTIMAL  KARHUNEN-LOEVE ) ' 

READ (  IN ,  3 )  TRCODE 
3  FORMAT(Il) 

DIM  *  M  +  1 
NCLASS  *  0 
PRECLS  =  0 

IF ( INCODE  .HE.  2)  GOTO  9 

READ ( QLQLUN , A )  TOTCLS 

6  READ ( GLDLUN , A  )  CLASS, SAMP 
IF (CLASS  .EU.  0)  GOTO  65 
DO  7  I  =  1,10 

READ ( OLDLUN, A )  ( L IB ( CLASS , SAMP , I , J ) , J=1 ,  1  +  1 ) 

7  CONTINUE 
GOTO  6 

9  NSAMP  =  0 
GOTO  30 

10  IF (  INCODE  .EQ.  2)  GOTO  65 

UR  ITE ( OUT , 15 )  NCLASS, NSAMP 

15  FORMAK IX, 'CLASS:  ',12,'  NUMBER  OF  SAMPLES:  ',12) 

UR ITE( OUT , 20 ) 

20  F0RMAK1X, 'FINISHED  UITH  SAMPLING?  (TYPE  Y  FOR  YES)  ',*) 

READ (IN, 25)  INST 
25  FORMAT ( A1 ) 

IE ( INST  .ECU  'Y'>  GOTO  100 

CALL  UIEU 

30  UR ITE ( OUT, 35) 

35  FORMAK IX, 'ENTER  CLASS  NUMBER:  ',*> 

RE AD (  IN, A)  CLASS 

IF ( CLASS  .EQ.  PRECLS)  GOTO  40 
NCLASS  *  NCLASS  +  1 
PRECLS  =  CLASS 
NSAMP  =•  0 
C 

C  CALCULATE  THE  AK  PARAMTERS  OF  THE  EDGES  OE  THE  TRAINING  OBJECTS 

C 

40  CALL  SNAP 

CALL  SCAN(BNDCI, BOUND, XPIX,YPIX) 

IF ( BNDCT  .EQ.  0)  GOTO  99 
DO  45  N  =  1, BNDCT 

T1  =  SECNDS (0.0) 

CALL  GETRAD(XPIX( 1 ,N> , YP IX ( 1 , N ), N , BOUND , RAD ) 

CALL  PARAM(N, RAD, THETA, ALL) 

DELTA  >  SECNDS ( T1 > 

UR  ITE ( DEV , A )  'TIME  FOR  GETRAD  AND  PARAM  IS:  '.DELTA 

45  CONTINUE 

C 

C  FORM  THE  SAMPLE  FEATURE  VECTORS  FROM  THE  AR  PARAMETERS 

C 

DO  60  J  ■  1, BNDCT 

NSAMP  ■  NSAMP  +  1 
DO  50  I  ■  1 , D IM 

X< I, NSAMP, CLASS)  ■  THETA ( I , J ) 


CONTINUE 

CONTINUE 

CALCULATE  IHE  ORIENTATION  OF  THE  EDGES 

CALL  ROTATE(BNDCT, BOUND, ORIENT) 

CALL  OUTPUT (8NDCT, BOUND, RAD, ORIENT) 

GOTO  99 

UR ITE ( OUT, 67 ) 

FORMAT( IX, 'FINISHED  UITH  TRAINING?  (TYPE  Y  FOR  YES):  '*> 
READ <  IN , 68 )  INST 
FORMAT ( A1 ) 

IF (  INST  .EQ.  'Y'>  GOTO  100 
UR ITE  <  OUT, 70  ) 

FORMATUX,  '  INPUT  CLASS  NUMBER:  '*> 

READ (  IN , SO )  CLASS 
UR ITE ( OUT , 75 ) 

FORMAT( IX, ' INPUT  DESIRED  NUMBER  OF  SAMPLES:  '$) 

READ (  IN , 90 )  NSAMP 
FORMAK  13) 

DO  90  SAMP  =  1, NSAMP 
DO  85  I  =  1 , DIM 

X( I, SAMP, CLASS)  =  LIB(CLASS,SAMP,M, I) 

CONTINUE 

CONTINUE 

IE (CLASS  .EQ.  PRECLS)  GOTO  99 
NCLASS  =  NCLASS  +  I 
PRECLS  =  CLASS 

LOOP  BACK  TO  GET  MORE  SAMPLES 
GOTO  10 

TRANSFORM  THE  FEATURE  VECTORS 


GOTO (110, 130, 130)  TRCODE 

CALL  MESELtX, NS AMP, NCLASS, TRANS, T) 

GOTO  150 

CALL  KLSEL(X, NS AMP, NCLASS, TRANS, I) 

GOTO  150 

CALL  PREKL(X, NS AMP, NCLASS, TRANS,!) 

FIND  THE  AVERAGE  OF  THE  TRANSFORMED  SAMPLES 

DO  180  K  =  1, NCLASS 
DO  170  J  =  1 , D IM 
SUM  =0.0 

DO  160  I  =  1, NSAMP 

SUM  =  SUM  +  T(  I, J,K) 

CONTINUE 

TA(J,K>  =  SUM  /  FLOAT ( NSAMP ) 

CONTINUE 

CONTINUE 

URITE  THE  TRANSFORMED  SAMPLE  FEATURES  TO  A  TRAIN.DAT  FILE 


200 


URITE!NEULUN,A>  Dili 


DO  210  I  =  l.DIM 

US  ITECNEULUN ,  205 )  ( TRANS (  I , J > , J= 1 , D IM > 

203  F0RMAT!1X,<DIM>E11.5> 

210  CONTINUE 

UR  ITE (NEULUN,230 )  NSAMP, NCLASS 
230  FORMAT! IX, 213) 

DO  260  K  =  1 , NCLASS 

URITE(NEULUN,255>  ! TA C J , K ) , J  =  1 , D  IM > 

255  FORMAT  ( 1 X , D  IM>E  1 1.5) 

260  CONTINUE 

IFCDIM-2  .GT.  0)  GOTO  290 
UR  ITE ( DEO  ,235  >  (I, 1  =  1, DIM) 

285  E0RMAT< 'O' , 'THE  AK  PARAMETERS  OF  THE  TRAINING  DATA'//, 

A  '  CLASS  SAMPLE  ' , 7X , <D IM> < ' M ' ,  12 , 7X ) / ) 

GOTO  305 

290  UR ITE ( DEV , 300 )  (I, 1=1, DIM) 

300  FORMAT( '0' , 'THE  AR  PARAMETERS  OF  THE  TRAINING  DATA'//, 

A  '  CLASS  SAMPLE  ' , 7X , CD IM- 1> ! ' M ' ,  12 , 7X ) , ' M ' , 12 , 3X , / ) 

305  DO  350  K  =  1, NCLASS 

DO  340  3  =  1 ,NSAMP 

UR ITE ! DEV , 330  )  K,J,(X(I,J,K),I*1,DIM) 

330  FORMAT! IX, 2 < 13, 5X ) , <D  IM>F10 .5 > 

340  CONTINUE 

350  CONTINUE 

UR  ITE ( DEU , 400 )  (I, 1=1, DIM) 

400  FORMAT! '0' , 'THE  TRANSFORMED  PARAMETERS  OF  THE  TRAINING  DATA'// 

A  '  CLASS  SAMPLE  ' ,7X,<DIM>( ' T ' ,  12 ,7X  >/ ) 

DO  450  K  =  1, NCLASS 

DO  440  I  =  1 ,NSAMP 

UR ITE! DEV, 430 )K,  I , ( T (  I , J , K ) , J  =  1 , D IM ) 

430  FORMAT! IX, 2! 13 , SX ) , <D IM>F10 . 4 ) 

440  CONTINUE 

450  CONTINUE 

CLOSE ! UN  11  =  NEULUN) 

CLOSE! UN  IT=OLDLUN) 
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C A******* A*** ******** A* A A*** **************************** A*A***AA*A**AA 
C 
C 

c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 

c ********************************************************************* 


FILE  UORK1.FCIR  LAST  REVISION  JULY  1,1985 

ROUTINES  CONTAINED  IN  THIS  FILE  ARE: 

SETUP 

THRESH 

SCAN 

GEICEN 

INTPLT 

OBPLOT 

FINPLT 

RADPLT 


SETUP 


SUBROUTINE  SETUP 
C 

C  THIS  SUBROUTINE  INITIALIZES  THE  IMAGING  SYSTEM  HARDWARE  AND 

C  PROMPTS  THE  USER  FOR  DESIRED  VALUES 

C 

INTEGER  DEV, IN, OUT, STEP, INRAD, M 
COMMON  DEV, IN, OUT, STEP, INRAD, M 


CALL  SELGRP ( 1 > 
CALL  AP IN  IT 
CALL  FB IN  IT 
CALL  LU IN  IT 
CALL  SYNC  <  0 ) 
CALL  VIEW 


UR ITE ( OUT , 40  > 

40  FORMATdX,  '  INPUT  DESIRED  OUTPUT  DEVICE:  ',*) 

READ (  IN , A )  DEV 

UR ITE ( OUT , 50 ) 

50  FORMATdX,  '  INPUT  DESIRED  SCAN  STEP:  ',*) 

READ ( IN , A )  STEP 

UR ITE ( OUT ,100) 

100  FORMATdX, 'INPUT  NUMBER  OF  RADIUS  VECTORS:  ',*) 
READ (  IN , *  )  INRAD 


UR ITE  <  OUT , 110) 

110  FORMAT< IX, ' INPUT  THE  DESIRED  AR  MODEL  ORDER  NUMBER:  ',*) 
READ ( IN , A )  M 


RETURN 

END 


C - 

C  THRESH 

C - 


SUBROUTINE  THRESH 


noon  onnoon  o  n  o 


C  THTS  SUBROUTINE  SETS  THE  IMAGE  THRESHOLD  AND  FILLS  THE  FRAME  BUFFER 

C 

INTEGER  DEM, IN,0UI,ANS, YES , THRS , CHIHRS ,  I 
COMMON  DEV, IN, OUT 

THRS  »  120 
INTLO  =  0 
INIHI  =  255 

CHTHRS  *  0 

100  UR ITE<  OUT , 110)  THRS 

110  FORMAK IX, 'THRESHOLD  =  ',14,'  ADD:  ',*> 

READ (  IN , 120  )  CHTHRS 
120  FORMAK  13) 

IF ( CHTHkS  .£0.  0)  GOTO  200 
THRS  =  THRS  +  CHTHRS 
GOTO  100 

200  CALL  SELLUT (0,0) 

DO  300  I  =  0,255 

IF (  I  .LT.  THRU)  CALL  SETLUT ( I , INTLO ) 

IF  < I  .GE.  THRS)  CALL  3ETLUT (  I , INTH  I ) 

300  CONTINUE 

RETURN 

END 


SCAN 


SUBROUTINE  SCAN(BNDCT,BOUND,X,Y) 

THIS  SUBROUTINE  SCANS  THE  FRAME  BUFFER  FOR  ALL  OF  THE  THRESHOLDED 
EDGES  OF  THE  OBJECTS  IN  THE  IMAGE.  ONCE  AN  EDGE  IS  FOUND  IT  CALLS 
GETCEN  TO  CALCULATE  THE  CENTROID,  AND  GETRAD  TO  CALCULATE  THE  RADIUS 
VECTOR  LENGTHS. 


INTEGER  RAPIX,RPIXEL,XL,ZP,XN,YN,ZN 
INTEGER  BNDCT,BOUND< 10,10) 

INTEGER  X< 2500, 10) ,Y(2500,10> ,Z,NUMPIX 
INTEGER  DEV, IN, OUT, STEP, INRAD 
INTEGER  FIRSTX,LASTX, INCX,STOPX 
INTEGER  F  IRSTY , L  ASTY 
INTEGER  EDGE, BCKGND, OBJ, XIEST 
INTEGER  BUF(0:S11) 

CHARACTER  INST 

LOGICAL  YESPLT , GOTONE 

COMMON  DEV, IN, OUT, SIEP,  INRAD 

FOR  WHITE  OBJECT  ON  BLACK  BACKGROUND:  BCKGND  »  0  AND  OBJ  =  255 
FOR  BLACK  OBJECT  ON  WHITE  BACKGROUND:  BCKGND  »  255  AND  OBJ  »  0 

BCKGND  =  0 
OBJ  »  255 

EDGE  =  255 
FIRSTX  =  0 


non  ono  onoonn  o  n  r>  ooon 


LASTX  =  511 
INCX  =  STEP 

STOPX  =  (LASTX- FIRSTX+1)  /  STEP 
FIRST*  »  0 
LASTY  =  479 
BNDCT  «  0 

UR ITE ( 5 , 50 ) 

50  FORMAT ( IX , ' TYPE  'Y1  IE  THE  OBJECT  PLOT  IS  DESIRED:  ',*> 

READ  ( 5,55)  INST 
55  FORMAT  <  A1 ) 

YESPLI  =  .FALSE. 

IF ( INST  .NE.  'Y')  GOTO  100 
YESPLT  *  .TRUE. 

CALL  INTPLT 

100  T1  =  SECNDS (0.0) 

DELTA  =  0.0 

DO  500  YN  =  FIRSTY, LASTY, STEP 

READ  EVERY  INCX  FROM  YN  LINE  IN  THE  FRAME  BUFFER  INTO  BUFF 

CALL  RSCAN(YN,FIRSTX, LASTX, INCX,BUF> 

ZP  -  BUF(O) 

DO  400  XL  *  1, STOPX 
ZN  -  BUF ( XL ) 

CHECK  FOR  AN  EDGE  AND  TRACE  BACK  FROM  HOP  TO  REAL  EDGE 

FOR  WHITE  ON  BLACK  USE  (ZN-ZP)  IN  NEXT  STATEMENT 
FOR  BLACK  ON  WHITE  USE  (ZP-ZN)  IN  NEXT  STATEMENT 


IF( (ZN-ZP)  .NE.  EDGE)  GOTO  300 
XN  *  (STEP  *  XL)  ♦  FIRSTX 
150  XTEST  *  XN  -  1 

Z  -  RPIXEL(X1ESI,YN> 

IF ( Z  .GT.  BCKGND  .AND.  Z  .LE.  BNDCT)  GOTO  300 
IF ( Z  .EQ.  BCKGND)  GOTO  200 
XN  =  XTEST 
GOTO  150 

TURTLE  AROUND  THE  EDGE  AND  CALCULATE  THE  CENTROID 

200  N  *  BNDCT  ♦  1 

CALL  GETCEN(XN, YN,X( 1 ,N> ,Y(1,N) , BNDCT , BOUND , GOTONE ) 
IF ( . NOT .  GOTONE)  GOTO  300 

PLOT  OUT  EDGES 

IF ( . NOT .  YESPLT)  GOTO  300 
DELIA  »  DELTA  ♦  SECNDS(Tl) 

NUMPIX  ■  B0UND<6, BNDCT) 

CALL  OBPLOT(X,Y, NUMPIX) 

T 1  *  SECNDS(O.O) 

300  ZP  *  ZN 

400  CONTINUE 

500  CONTINUE 


noo  onooooo 


106 


DELTA  =  DELTA  ♦  SECNDS(Tl) 

C  IF ( YESPLT )  CALL  FINPLI 

WRITE (OUT, GOO )  DELTA 

600  FORMATUX, 'THE  SCAN  AND  TURTLE  TIME  IS:  'fF18.6f'  SEC'/) 

1000  RETURN 
END 


GETCEN 


SUBROUTINE  GETCEN < XS , YS , X , Y , BND .1 , BOUND , GOTONE ) 

THIS  SUBROUTINE  CALULATES  THE  CENTROID  OF  THE  BOUNDARY .  THE  EDGE 
IS  TRACED  USING  A  TURTLE  ALGORITHM. 

INTEGER  XS, YS , BNDCT , BOUND < 10, 10) 

INTEGER  XIEST 

INTEGER  RPIXEL,RAPIX,XT, YT,Z,BNDINT 
INTEGER  XREG,YR£6, NUMP IX , XC , YC 
INTEGER  X ( 2500 ) , Y ( 2500 ) 

INTEGER  N , TOL , BCKGND 

LOGICAL  OBJECT, RIGHT, UP, GOTONE 


TOL  *  25 
BCKGND  ■  0 
N  *  BNDCT  +  1 
BND INT  *  N 

CALL  UPIXEL(XS,YS, BND INT ) 

XREG  =  XS 
YREG  =  YS 
NUMP IX  =  1 
X( NUMP IX)  =  XS 

Y  <  NUMP IX )  =  YS 
XT  =  XS  -  1 

YT  =  YS  -  1 
UP  =  .TRUE. 

ISO  Z  =  RP IXEL ( XT , YT ) 

IF  ( Z  .EO.  BCKGND)  GOTO  250 

IF( (XT.EQ.XS)  .AND.  (YT.EO.YS)  .AND.  ( NUMP IX . GT . 10 ) )  GOTO  500 
IF ( Z  .LI.  BND INT )  GOTO  500 
200  OBJECT  =  .TRUE. 

IF< (XT.EQ.X(NUMPIX) ) .AND. (YI.EQ.Y (NUMPIX) ) >  GOTO  260 
XREG  =  XREG  ♦  XT 
YREG  =  YREG  +  YT 
NUMPIX  =  NUMPIX  ♦  1 
X ( NUMPIX )  -  XT 

Y  <  NUMP IX )  =  YT 

CALL  UPIXEL(XT,YT,BNDINT) 

GOTO  260 

250  OBJECT  «  .FALSE. 

MOVE  UP  OR  DOWN 
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260 


?7e 


290 


300 


350 

360 

375 

390 

500 


IF< (UP. AND. OBJECT) .OR. ( ( .NOT. UP) .AND. ( .NOT. OBJECT) ) )  GOTO  275 
XT  =  XT  +  1 
RIGHT  =  .TRUE. 

GOTO  290 
XT  »  XT  -  1 
RIGHT  *  .FALSE. 

CHECK  FOR  OBJECT  PIXEL 

Z  =  RP IXEL  ( XT , YT ) 

IF  (  Z  .£0.  BCKGND)  GOTO  350 

IF < <XT.EQ.XS)  .AND.  (YT.EQ.YS)  . AND . ( NUMP IX . GT . 1 0 ) )  GOTO  500 
IF ( Z  .LI.  BNOINT)  GOTO  500 
OBJECT  =  .TRUE. 

IF <  ( XT . EQ . X ( NUMP IX ) )  .AND.  < YT . EO . Y ( NUMP IX > > >  GOIO  360 

XREG  =  XREG  +  XT 

YREG  =  YREG  +  YT 

NUHP  IX  =  NUMP  IX  +  1 

X ( NUMP  IX )  =  XT 

Y ( NUMP  IX )  =  YT 

CALL  UPIXEL(XT,YT,BNDINT) 

GOTO  360 

OBJECT  =  .FALSE. 


OR. ( ( . NOT. RIGHT) .AND. ( .NOT. OBJECT) )  ) 


MOVE  UP  OR  DOWN 

IF  < (RIGHT. AND. OBJECT) . 
GOTO  375 
YT  «  YT  ♦  1 
UP  *  .FALSE. 

GOTO  390 
YT  =  YT  -  1 
UP  =  .IRUE. 

GOTO  190 


CALCULATE  THE  XC,YC  OF  THE  CENTROID  PIXEL 

IF ( NUMP IX  .GT.  TOL)  GOTO  600 
GOTONE  =  .FALSE. 

GOTO  1000 


600  BNDCT  =  BNDCT  +  1 
GOTONE  =  .TRUE. 

XC  =  XREG  /  NUMP IX 
YC  =  YREG  /  NUMP  IX 

X ( NUMP IX+ 1 )  =  XS 
Y  <  NUMP  IX+ 1 >  =  YS 
BOUND ( 1 , N )  =  BNDINT 
BOUND ( 2 ,  N )  =  XS 
BOUND ( 3 , N )  *  YS 
BOUND ( 4 , N )  =  XC 
BOUND <  5 ,  N  >  =  YC 
BOUND  <  6 ,  N )  =  NUMP IX 


1000 


RETURN 

END 


INTPLT 


c 


SUBROUTINE  INIPLI 
C 

C  THIS  SUBROUTINE  INITIALIZES  THE  PLOTTER. 

C 

INIEGER  LUN 
COMMON  /DEV  DAT/LUN 

LUN  =  5 
CALL  IXTCLR 
CALL  PLTINT 
CALL  PLTCLR 
C 

C  PLOT  THE  FRAME 

C 

CALL  PLTCUR (0,0) 

CALL  PLTVEC (511,0) 

CALL  PLTUEC ( 511, A79  > 

CALL  PLTVEC (0,479) 

CALL  PLTVEC(0,0> 

RETURN 

END 


C - 

C  OBPLOT 

C - 

SUBROUTINE  OBPLOT ( X , Y ,NUMP IX ) 

C 

C  THIS  SUBROUTINE  PLOTS  OUT  THE  POINTS  COMPRISING  THE  EDGE  COORDIN- 

C  ATES  HELD  IN  ARRAY  X  AND  ARRAY  Y. 

C 

INTEGER  X<3000> , Y(3000> ,NUMPIX 


DO  100  I  =  1 , NUMP  IX 
IX  =  X( I) 

IY  =  Y<  I) 

CALL  PLTPNT ( IX,  IY) 
100  CONTINUE 

RETURN 

END 


C - 

C  FINPLT 

c - 

SUBROUTINE  FINPLT 
C 

C  THIS  SUBROUTINE  PRINTS  OUT  A  HARD  COPY  OF  THE  PLOT  IF  THE  USER  SO 

C  DESIRES  IT. 

C 

CHARACTER  INST 
UR ITE ( 5 , SO ) 


wl 


50  FORMAT ( IX , 'TYPE  *Y’  IF  HARD  COPY  IS  DESIRED:  ',$> 

READ ( 5,55)  INST 
55  FORMAT  <  A1 ) 

IF ( INST  .NE.  'Y'>  GOTO  100 
CALL  HARDCOPY 
100  CALL  PLTCLR 

CALL  TXTCLR 
CALL  PLTOFF 

RETURN 

END 

c - 

C  RADPLT 

c - 


SUBROUTINE  R ADPLT ( RAD , EDGE > 

INTEGER  EDGE , NUMRAD 
REAL  RAD ( -2 : 250 ,  10 ) 

N  =  EDGE 

MAXY  =  479 

RNUMRAD  =  RAD(-2,N> 

NUMRAD  =  INT ( RAD ( -2 , N  >  ) 

INC  =  512  /  NUMRAD 
IXB  =  0 
IYB  =  MAXY 

DO  100  I  *  1 , NUMRAD 

IXB  =  INT ( FLOAT ( I)A511. O/RNUM RAD) 
IYT  =  MAXY  -  ( INT ( RAD ( I ,N  >  >  A  2> 
CALL  PLTCUR < IXB, IYT) 

CALL  PLTVEC( IXB, IYB) 

100  CONTINUE 

RETURN 

END 


no 


CAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAftA 

C 

C  FILE  RADI. FOR  LAST  REVISION  JULY  I,  1985 

C 

C  THIS  FILE  CONTAINS  THE  SUBROUTINE: 

C 

C  GETRAD 

C 

CAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA 

c - 

C  GETRAD 

c - 

SUBROUTINE  GETS AD ( X , Y , EDGE , BOUND , RAD ) 

C  THIS  SUBROUTINE  FINDS  AND  CALCULATES  THE  RADIUS  VECTORS  OF 

C  AN  OBJECT  EDGE.  IT  USES  AN  EQUI-ANGLE  ALGORITHM  SIMILAR 

C  TO  DUBOIS  I1984J.  THIS  SUBROUTINE  READS  THE  EDGE  PIXEL  LOCATIONS 

C  WHICH  ARE  TEMPORARILY  STORED  IN  THE  X  AND  Y  ARRAYS. 

C 

INTEGER  XC2500) , Y ( 0500) , BOUND C 10,10), EDGE 
INTEGER  DEV, IN, OUT, I , N , M , XC  ,  YC , NUMP  IX 
INTEGER  INRAD, NSLOPE,LOUVEC,HIVEC, SEC 
INTEGER  XMAX , Y MAX 

REAL  RAD< -2:250,10) , SLOPE (0:90) ,P I, PH  I 

REAL  VEC,XDIF, YD  IF, YTRY1, YTRY2,YDR1 ,YDR2 

REAL  YD  1 , YD2 , MAXRAD 

LOGICAL  START, STOP 

COMMON  DEV, IN, OUT, STEP, INRAD 

DATA  PI/3.141592654/ 

N  *  EDGE 

XC  =  BOUND ( 4 , N ) 

YC  =  BOUND ( 5 , N ) 

NUMP  IX  =  BOUND ( 6 , N ) 

NSLOPE  =  INRAD/4  -  1 

PHI  =  <2.0  A  PI)  /  FLOAT ( INRAD ) 

DO  50  I  =  0, NSLOPE 

SLOPE  < I )  =  T  AN ( PH  I  A  FLOAT (  I ) )  A  (5. 0/4.0) 

50  CONTINUE 

SLOPE ( NSLOPE  ♦  1)  =  1250.0 

STOP  =  .FALSE. 

START  =  .TRUE. 

I  =  0 

PREVEC  =  0 
SEC  =  0 
M  -  0 

MAXRAD  =  0.0 
C 

C  CALCULATE  THE  SLOPE  OF  THE  LINE  FROM  THE  EDGE  PIXEL  TO  THE  CENTROID 
C 

100  1=1*1 

IF (  I  .GE.  NUMP IX  >  STOP  =  .TRUE. 

XD IF  =  FLOAT ( X (  I )  -  XC ) 

YD  IF  =  FLOAT ( Y (  I )  -  YC) 

IF ( XD  IF  ,EU.  0.0)  GOTO  700 
VEC  =  ABS ( YD IF/XD  IF ) 

IF ( ABS ( VEC  -  PREVEC)  .LE.  .0001)  GOTO  100 


•  '  -.V/  O  O  %■  v* 


LOUVEC  =  SEC 
PREVEC  =  VEC 
C 

C  FIND  OUT  WHICH  SECTOR  WE  ARE  IN 
C 

110  IE ( VEC  .GE.  SLOPE ( LOUVEC )  .AND.  VEC  .LE.  SLOPE ( LOW VEC + 1 ) ) 
k  GOTO  120 

IF ( VEC  .LE.  SLOPE ( LGWVEC  i  .AND.  VEC  .GE.  SLOPE ( LuWVEC-1  )  ) 
k  GOTO  140 

IF  <  LOWVEC  .GE.  NSLOPE)  GOTO  115 
LOWVEC  =  LOUVEC  +  1 
GOTO  110 

.115  LOUVEC  =  1 

GOTO  110 

ICO  H I VEC  =  LOUVEC  +  1 

GOTO  200 

140  HIVEC  =  LOUVEC  -  1 

C 

C  CALCULATE  THE  REFERENCE  POSITION  BETWEEN  LOVEC  AND  HIVEC 
C 

200  YTRY1  =  ABS(XDIF)  *  SLOPE ( LOUVEC ) 

YTRY2  =  ABS(XDIF)  *  SLOPECHIVEC) 

YDR1  =  ABS  <  YD  IF  >  -  YTRY1 
YDR2  =  ABS  <  YD  IF )  -  YTRY2 
C 

C  LOOK  AT  THE  NEXT  POSITION  ALONG  THE  EDGE 
C 

300  1=1+1 

IF ( I  .GE.  NUMPIX)  SIOP  =  .TRUE. 

XD IF  =  FLOAT ( X (  I)  -  XC) 

YD  IF  =  FLOAT ( Y ( I)  -YC) 

IF ( XD  IF  .EQ.  0.0)  GOTO  700 
GOTO  800 

700  HIVEC  =. NSLOPE  +  1 

LOUVEC  =  NSLOPE 
GOTO  1000 
C 

C  CALCULATE  THE  NEU  POSITION  WITH  RESPECT  TO  HIVEC  AND  LOUVEC 
C 

800  YTRY1  =  ABS(XDIF)  *  SLOPE < LOUVEC ) 

YTRY2  =  ABS(XDIF)  *  SLOPE(HIVEC) 

YD  1  =  ABS ( YD  IF  >  -  YTRY1 
YD2  *  ABS(YDIF)  -  YTRY2 
C 

C  CHECK  IF  UE  HAVE  CROSSED  LOVEC 
C 

IF ( ( YD1  .GE.  0.0  .AND.  YDR1  .LE.  0.0)  .OR.  ( YD1  .LE.  0.0 
k  .AND.  YDR 1  .GE.  0.0>>  GOTO  800 
C 

C  CHECK  IF  UE  HAVE  CROSSED  HIVEC 
C 

IF ( ( YD2  .GE.  0.0  .AND.  YDR2  .LE.  0.0)  .OR.  ( YD2  .LE.  0.0 
k  .AND.  YDR2  .GE.  0.0>>  GOTO  1000 
C 

C  HAVEN'T  CROSSED  LOVEC  OR  HIVEC  YET,  SO  GET  THE  NEXT  EDGE  PIXEL 
C 


IF ( STOP )  GOTO  2000 
GOTO  300 
C 

C  YEH,  MADE  IT  ACROSS  A  SLOPE  MARKER  CLOVEC  OR  HIVEC) 
C  NOU  UE  CAN  CALCULATE  THE  RADIUS  VECTOR  LENGTH 
C 

C  UE  CROSSED  LOVEC 

C 

900  IF ( LOUVEC  .EG.  SEC)  GOTO  100 

M  =  M  +  1 

RAIKM.N)  =  5  Q  R  T  ( (XDIFA1.2j)A*2  ♦  YDIFAA2) 

SEC  =  LOUVEC 
IF  (  STOP  )•  GOTO  2000 
GOTO  1200 
C 

C  UE  CROSSED  HIVEC 

C 

1000  IF ( H IVEC  .EG.  SEC)  GOTO  100 

M  =  M  ♦  1 

R AD ( M , N )  =  SGRT( (XDIF*1.25)**2  +  YDIF*A2) 

SEC  =  HIVEC 
IF  (STOP)  GOTO  2000 

1200  IF(RAD(M,N)  .LE.  MAXRAD )  GOTO  100 
MAXRAD  =  R AD  <  M , N ) 

XMAX  =  X (  I ) 

YMAX  =  Y<  I) 

GOTO  100 

2000  RAD<-2,N)  =  FLOAT ( M ) 

RAD(-1,N)  =  MAXRAD 
SOUND  <  9 ,  N )  =  XMAX 
BOUND( 10, N)  =  YMAX 
C  CALL  FINPLT 

C  CALL  INTPLT 

C  CALL  RADPLT(RAD.EDGE) 

C  CALL  FINPLT 

RETURN 

END 


CAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA 

C 

C  FILE  RAD2.F0R  LAST  REVISION  JULY  26,1985 

C 

C  SUBROUTINES  CONTAINED  IN  THE  FILE: 

C 

C  GETRAD 

C 

CAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA 

SUBROUTINE  GETR AD < X , Y , BNDCT , BOUND , RAD ) 

C 

C  THIS  SUBROUTINE  CALCULATES  THE  OBJECT  EDGE  RADIUS  LENGTHS 

C  USING  THE  EQUAL  ARC  LENGTH  ALGORITHM  . 

C 

INTEGER  DEV, IN, OUT, STEP, INRAD 

INTEGER  X ( 2500 > , Y ( 2500 ) , BNDCT , BOUND ( 1 0 , 1 0 > 

INTEGER  XC, YC,NUMPIX, INC 
INTEGER  XMAX , YMAX 

REAL  RAD < -2:300, 10) ,XD  IF,  YD  IF, MAXRAD 

COMMON  DEV, IN, OUT, STEP, INRAD 

N  =  BNDCT 

XC  *  BOUND ( 4 ,N) 

YC  =  BOUND ( 5, N> 

NUMPIX  =  BOUND<  6,N) 

MAXRAD  =  0.0 

DO  200  I  *  1,  INRAD 

INC  =  I  A  NUMPIX  /  INRAD 
XD  IF  =  FLOAT ( X (  INC )  -  XC) 

YD  IF  *  FLOAT ( Y < INC )  -  YC) 

RAD ( I , N )  =  SQRT( (XDIFA1 .25>AA2  ♦  YDIFAA2) 

CALL  PLICUR< XC , YC ) 

CALL  PLTVEC(X( INC) , Y( INC) ) 

IF ( R AD  < I , N )  .LI.  MAXRAD)  GOTO  200 
MAXRAD  =  RAD(  I,N) 

XMAX  =  X<  INC) 

YMAX  =  Y ( INC ) 

CONTINUE 

RAD  < -2 , N )  =  FLOAT(INRAD) 

R AD  < - 1 ,  N )  =  MAXRAD 
BOUND  <  9 , N  )  =  XMAX 
BOUND ( 10 , N )  =  YMAX 
C  CALL  FINPLT 

C  CALL  INTPLI 

C  CALL  RADPLT (RAD , N > 

C  CALL  FINPLT 

RETURN 
END 


C 

C 


200 


CAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA 

C 

C  FILE  PARAM . FOR  LAST  REVISION  JULT  1,1985 

C 

C  THIS  FILE  CONIAINS: 

C 

C  TEST  ROUTINE  FOR  PARAM 

C  PARAM 

C 

caaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaa 

c 

C  THIS  IS  A  TEST  ROUTINE  USED  TO  CHECK  THE  AR  COEFF- 

C  ICIENTS  AND  THE  SUBROUTINE  PARAM  WHICH  CALCULATES 

C  THE  AUTOCORRELATION  AND  THE  AUTOREGRESSIVE  COEFFICIENTS 

C  OF  THE  TIME  SERIES. 

C 

C  THE  TEST  ROUTINE 


INTEGER  DEV 

REAL  RAD <-21300, 10) ,ALPH( 11 , 10) ,RX(0: 10, 10) 

COMMON  DEV 
DEV  =  5 


N  =  1 

R AD ( 1 , N )  =  1. 

RAO  <  2 , N )  =  1.44 

RAO  <  3 , N )  =  1 . 44ARA0 ( 2 , N  >  -  1.21 
DO  10  J  «  4,40 

RADU.N)  *  l.44*RAUU-l,N) 
CONTINUE 


1.2&ARADC J-2,N)  ♦  .  8 1  ARAD  (  3-3  ,  N  ) 


3 

|  ® 
•Ki 


RAD<-2,N)  *  40. 

CALL  PARAM < N , RAD, ALPH , RX ) 


SUBROUTINE  PAR AM < EDGE , RAD , THETA , ALL ) 

THIS  SUBROUTINE  CALCULATES  THE  FEATURE  VECTOR  PARAMETERS  OF  THE 
EDGES  OF  THE  OBJECT.  THE  PARAMETERS  ARE  THE  AR  COEFFICIENTS  AND 
THE  TERM  ALPHA/SORT < BETA ) . 

INTEGER  DEV, IN, OUT, STEP, INRAD, M 
INTEGER  EDGE , NUMRAD 

REAL  SUM , SUM1 , SUM2 

REAL  ALPHA, BETA, ALBET 

REAL  RAD <-2:250, 10) , THET A < 1 1 , 10 ) , RX < 0 : 10 , 1 0 > 

REAL  ALL< 10, 11 , 10) 

REAL  RXN<0:10,10> 

REAL  SIGMASC10) ,A<10,10) 

COMMON  DEV, IN, OUT, STEP, INRAD, M 


N  s  EDGE 
DO  30  J  =  1,10 

DO  20  I  =  1,10 


M 

mmmmm&mm-  mm 

A<  1,3)  *  0.0 
20  CONTINUE 

30  CONTINUE 

NUMRAD  »  1NI<RAD<-2,N> ) 

C 

C  CALCULATE  THE  MEAN 

C 

SUM  *  0.0 

DO  SO  I  =  1, NUMRAD 

SUM  =  SUM  *  R  AD  < I,N> 

30  CONTINUE 

MEAN  =  SUM  /  FLOAT ( NUMRAD  > 

C 

C  CALCULATE  THE  AUTOCORRELATION  COEFFICIENTS 

C 

DO  100  I  =  0 ,  M 
SUM  =  0.0 

DO  733=  1, NUMRAD 
K  =  3  +  I 

IF ( K  .GT.  NUMRAD)  K  =  K  -  NUMRAD 

SUM  =  SUM  +  ( <  R AD ( 3 , N )  -  MEAN)  *  ( RAD  <  K , N )  -  MEAN)) 
75  CONTINUE 

RX<I,N>  =  SUM  /  FLOAT ( NUMRAD ) 

RXN  <  I  f  N  )  =  RX ( I ,  N  >  /  RX ( 0 , N ) 

100  CONTINUE 
C 

C  CALCULATE  THE  AUTOREGRESSIVE  COEFFICIENTS  USING  THE  DURBIN 

C  RECURSIVE  ALGORITHM 

C 

A ( 1 , 1 )  *  -RX<1,N>  /  RX ( 0 , N ) 

S  IGMAS  < 1 )  -  <1.0  -  ( A ( 1 r 1 >  A  A( 1 , 1 ) ) >  A  RX ( 0  p  N  > 

DO  300  K  =  2, M 
SUM  =  0.0 

DO  230  3  =  1 »K  -  1 

SUM  =  SUM  ♦  < A (K-l , 3)  A  RX(K-3,N)> 

250  CONTINUE 

A ( K , K )  =  -<RX<K,N)  +  SUM)  /  S  IGMAS < K-l) 

DO  275  I  =  1 f  K - 1 

A<K,I)  =  A < K- 1,1)  +  <A<K,K)  A  A<K-1,K-I>> 

273  CONTINUE 

S  IGM AS  <  K )  =  < 1 . 0- <  A ( K , K )  A  A ( K  t  K ) ) )  A  SlGMAS(K-l) 

300  CONTINUE 

C 

C  FORM  THE  FEATURE  VECTOR  FOR  THE  EDGE  FROM  THE  LAST  COLUMN 

C  OF  THE  MATRIX  A.  THIS  COLUMN  CORRESPONDS  TO  THE  DESIRED 

C  AR  MODEL  ORDER  NUMBER.  ALSO  FORM  THE  SUM  OF  THE  PARAMETERS 

C  NEEDED  IN  THE  CALCULATION  OF  ALPHA. 

C 

SUM  *  0.0 
DO  400  K  =  1 ,  M 

THETA ( K , N  >  =  A ( M , K  > 

SUM  =  SUM  +  THETA ( K , N ) 

400  CONTINUE 

ALPHA  »  MEAN  A  <1.0  -  SUM) 


oooo  noon  n  n  n 


c 

c 


410 

420 

C 

C 

C  425 
C 

C  450 

700 

702 
705 

710 
715 


CALCULATE  BETA 

BETA  =  0.0 
DO  420  I  =  1 , NUMRAD 
SUM  =  0.0 
DO  410  J  =  1,M 
K  =  I  +  J 

IF ( K  .GT.  NUMRAD)  K  =  K  -  NUMRAD 
SUM  =  SUM  +  THETA ( J ,  N )  A  R  A  D  (  K ,  N  > 
CONTINUE 

SUM1  =  RAD( I, N)  -  ALPHA  -  SUM 
BETA  =  BETA  +  SUM  1 A A2 
CONTINUE 

BETA  =  BETA  /  FLOAT < NUMRAD ) 

BETA  =  S IGMAS (  M  > 

ALBEI  =  ALPHA  /  SORT ( BETA ) 

THETA  <  M+ 1 , N )  =  ALBET 


UR  HE  <  DEV , 425 )  (I,(A(K,I),K=1,M),I=1,M) 
FORMAT  < IX , 13 , 3X , <M>  F8 . 3 ) 

URITE(DEU, 450)  < J , S IGMAS < J ) , J=1 , M ) 
FORMATUX,  I3.E18.6) 


UR ITE ( DEV , 700 ) 

FORMAT ( IX , / '  THE  AUTOCORRELATION  COEFFICIENTS  (PARAM):  '/) 

UR ITE ( DEV ,705 )  (J, ( RX ( J  ,  N ) fN  =  l,6NDCT) , J  =  0,M) 

UR  ITE ( DEV , 702 ) 

FORMAT ( IX , /  '  THE  NORMALIZED  AUTOCORRELATION  COEFS  (PARAM):  '/) 
UR ITE (DEV ,705)  (  J , ( RXN ( J , N ) ,N=1,BNDCT) ,J  =  0,M) 

FORMAT (IX, 13 , 3X , <BNDCT>F10 . 5 ) 


URITE(DEV,710) 

FORMAT (IX,/'  THE  AUTOREGRESSIVE  COEFFICIENTS  (PARAM):  '/) 

UR ITE  <  DEV , 715  )  ( J ,( THETA ( J , N ), N= 1 , BNDCT ), J  =  1,M+1> 

FORMAT ( IX, 13 , 3X , • BNDCT>F1 0 . 5 ) 

i 

j 

RETURN 

END 
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caaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaa 

c 

C  FILE  KLSEL. FOR  LAST  REVISION  JULY  4,  1985 

C 

C  THIS  FILE  CONTAINS  SUBROUTINE: 

C 

C  KLSEL 

C 

caaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaa 

c - 

C  KLSEL 

C - 

c 

SUBROUTINE  KLSEL ( X , NS AMP , NCLASS , TRANS , Y ) 

C 

C  THIS  SUBROUTINE  CALCULATES  THE  KARHUNEN-LOEVE  EXPANSION 

C  TRANSFORM  MATRIX  AND  TRANSFORMS  THE  SAMPLES  OF  EACH 

C  PATTERN  CLASS  OF  THE  TRAINING  SET. 

C 

INTEGER  DEV, IN, OUT, STEP,  I N R A H , M 
INTEGER  OIM,JOBN, 12,  IER 
INTEGER  NCLASS, NSAMP, MAXI 
REAL  A(IOO) ,UK( 100) ,D( 11 )  ,Z( 100) 

REAL  X(11,25,10>,XMNU1,25,10>,R(11,11,10>,C<11,11> 

REAL  MEAN<11,10),MMT(11,11,10) 

REAL  SCALE, MAXE 

REAL  TRANS< 11,11) ,Y<25, 11,10) 

COMMON  DEV, IN, OUT, STEP, INRAD, M 

DATA  R/1210A0.0/ 

DATA  C/121A0.0/ 

T1  =  SECNDS  <  0.00 ) 

DIM  =  M  +  1 

CALCULATE  THE  MEAN  SAMPLE  FOR  EACH  CLASS 

SCALE  =  1.0  /  FLOAT ( NSAMP ) 

DO  50  N  =  1, NCLASS 
DO  40  I  =  1 , D IM 
SUM  =  0.0 
DO  30  J  =  1, NSAMP 

SUM  =  SUM  +  XC  I, J,N> 

CONTINUE 

MEAN (  I,N>  =  SUM  A  SCALE 
CONTINUE 
CONTINUE 

SUBTRACT  IHE  MEAN  FROM  EACH  SAMPLE  IN  THE  CLASS 

DO  100  N  =  1, NCLASS 
DO  75  K  =  l, NSAMP 
DO  SO  I  =  1 ,DIM 

XMN ( I , K , N  )  =  X  < I , K , N )  -  MEAN ( I , N ) 

CONTINUE 
CONTINUE 
CONTINUE 


nooooooonn  non 


C  FORM  THE  COVARIANCE  MATRIX  OF  EACH  CLASS 

C 

DO  COO  N  *  1 , NCLASS 

DO  1 75  K  =  X , NSAMP 
DO  150  I  *  I, DIM 

DO  125  J  =  1 ,  B  IM 

R(I,J,N>  =  R(I,J,N>  ♦  XMN<I,K,N>  *  XMN<J,f(,N> 

125  CONTINUE 

150  CONTINUE 

175  CONTINUE 

200  CONTINUE 

C 

C  CALCULATE  THE  AVERAGE  COVARIANCE  MATRIX 

C 

SCALE  =  1.0  /  FLOAT (NCLASS  A  NSAMP) 

DO  300  I  =  1 , D  I M 

DO  275  J  =  l.DIM 

DO  250  N  »  1, NCLASS 

C(I,J>  =  C(I,J)  +  (SCALE  A  R ( I , J , N ) ) 

250  CONTINUE 

275  CONTINUE 

300  CONTINUE 

C  DO  325  I  =  1,DIM 

C  UR ITE( DEV , A )  'C  ' ,  I ,  (C( 1 , 3 ) , J  *  l.DIM) 

C  325  CONTINUE 

C 

C  PUT  THE  LOWER  IRAINGLE  OF  MATRIX  C  INTO  ARRAY  A  (SYMMETRIC  STORAGE 

C  MODE) 

C 

K  =  1 

DO  340  I  =  1  ,  D IM 
DO  335  J  =  1,1 

A  <  K  )  =  C  <  I ,  J  > 

C  UR ITE ( D  E V , A )  'A  '  f  A  <  K ) 

K  =  K  ♦  1 

335  CONTINUE 

340  CONTINUE 

CALCULATE  THE  EIGENVALUES  AND  THE  EIGENVECTORS 

JOHN  =  2 
IZ  =  DIM 

CALL  EIGRS(A,DIM,JOBN,D,Z, IZ.UK, IER) 

K  »  1 

DO  350  I  *  l.DIM 
UR  ITE ( DEV , A  ) 

UR  ITE ( DEV , A )  'EIGENVALUE:  D(I> 

UR  ITE ( DEV , A ) 

DO  348  N  =  l.DIM 

UR  ITE ( DEV . A  >  'EIGENVECTOR:  ',Z(K> 

K  ■  K  ♦  1 

348  CONTINUE 

350  CONTINUE 


C  UR ITE (  5 ,  A  )  'UK(l)  ' , UK  < I ) 

C 

C 

C  FORM  THE  COLUMNS  OF  THE  TRANSFORMATION  MATRIX  FROM  THE  EIGENVECTORS 

C  CORRESPONDING  TQ  THE  LARGEST  EIGENVALUES. 

C 

IMAX  =  DIM  A  <  D IM  -  1) 

DO  S25  J  =  l.DIM 

DO  SIS  I  *  I , D IM 

TRANS ( I,J>  =  Z< IMAX+I) 

515  CONTINUE 

IMAX  =  IMAX  -  DIM 
525  CONTINUE 

C  UR ITE ( DEV , A  >  'THE  TRANSFORMATION  MATRIX' 

C  UR ITE  <  DEV , 527 )  < < TRANS < I , J > , J  =  1 , D  IM ) ,  1  =  1 , D IM > 

C  527  FORMAT ( IX , <D IM>FI5. 6 ) 

C 

C  TRANSFORM  THE  ORIGINAL  SAMPLES  OF  EACH  CLASS 

C 

DO  600  K  =  l.NCLASS 

DO  575  J  =  1 , NSAMP 
DO  550  N  =  1 , D  IM 
Y(J,N,K>  =  0.0 
DO  530  I  =  1 , D  IM 

Y  (  3  ,  N ,  K  >  »  Y ( J  ,  N  ,  K  )  +  (TRANSIT, N)  A  X(I,J,K)> 

530  CONTINUE 

550  CONTINUE 

575  CONTINUE 

600  CONTINUE 

DELTA  »  SECNDS  <  T1 ) 

UR  ITE ( DEV, A )  ' KLSEL  EXECUTION  TIME  IS:  '.DELTA 
C  UR ITE ( DEV , A )  'THE  TRANSFORMED  SAMPLES' 

C  URITE(DEV,700)  ( ( ( Y ( J , N , K ) , N= 1 , D IM ) , J  =  I, NSAMP), K  =  l.NCLASS) 

C  700  FORMAT ( IX , <D IM>F15 .6 ) 

RETURN 
END 


CAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA 
C  MESEL 

caaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaa 

SUBROUTINE  MESEL  IS  EXACTLY  IHE  SAME  AS  SUBROUTINE  KLSEL  EXCEPT 
THAT  THE  ORDERING  QE  THE  EIGENVECTORS  IN  THE  RESULTING  TRANS¬ 
FORMATION  MAIRIX  IS  REVERSED.  TO  AVOID  REPETITION,  THE  LISTING 
HAS  BEEN  OMITTED. 

C - 

C - 
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CAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAwAAAAAAAAAAAA 

c 

C  FILE  PREKL. FOR  LAST  REVISION  JULY  &,1985 

C 

C  THIS  FILE  CONTAINS  SUBROUTINE: 

C 

C  PREKL 

C 

c 

caaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaa 

c - 

C  PREKL 

c - 

SUBROUTINE  PREKL ( X, NSAMP , NCLASS f TRANS , Y > 

INTEGER  DEV, IN, OUT, STEP, INRAD, M 
INTEGER  DIM, JOBN,  IZ,  IER 
INTEGER  NCLASS, NSAMP 
REAL  SCALE, SUM 

REAL  X< 1 1,25,20), TRANS (11,11),Y( 25, 11, 20) 

REAL  MEAN( 11 ,20! ,XMN( 11,25,20) , AVEMEAN ( 1 1 ) 

REAL  SU< 11 , 11 ) ,R< 11 ,11 ,20) 

REAL  A (100) ,Z< 121 > ,D< 11) ,UK< 150) 

REAL  8(11,11) 

REAL  SbSUB (11,11),S8PR1ME(11,11),SB(11,11) 

COMMON  DEV, IN, OUT, STEP, INRAD, M 

DATA  R/2420A0. 0/ 

DATA  SU/121A0.0/ 

DATA  SB/121  AO. 0/ 

T1  =  SECNDS ( 0 . 00 ) 

DIM  =  M  +  1 

SCALE  =  1.0  /  FLOAT  <  NSAMP  > 

C 

C  CALCULATE  THE  MEAN 

C 

DO  100  N  «  l, NCLASS 
DO  75  I  =  1,BIM 
SUM  =0.0 
DO  50  J  =  1, NSAMP 

SUM  =  SUM  +  X(  I, J,N) 

50  CONTINUE 

MEAN (  I , N  )  =  SUM  A  SCALE 
75  CONTINUE 

100  CONTINUE 

C 

C  CALCULATE  THE  MEAN  VECTOR  A  MEAN  VECTOR  (TRANSPOSED)  FOR  EACH 

C  CLASS 

C 

DO  200  N  =  1, NCLASS 

DO  175  K  =  1, NSAMP 
DO  150  I  =  l.DIM 

XMN(I,K,N>  »  X (  I , K , N )  -  MEAN (  I , N ) 

150  CONTINUE 

175  CONTINUE 

200  CONTINUE 


ooo  n  no  o  o  o  ooonono  non 


c 

c 


CALCULATE  PART  OF  THE  COVARIANCE  MATRIX  FOR  EACH  CLASS 


DO  300  N  =  1 , NCL ASS 

DO  275  K  =  1 , NS AMP 
DO  250  I  =  1 , D IM 

DO  225  J  =  1 , D  IM 

R ( I ,  J  ,  N  )  =  R ( I ,  J  ,  N  )  +  XMN(I,K,N)  A  XMN (  J  ,  K  ,  N  ) 
225  CONTINUE 

250  CONTINUE 

275  CONTINUE 

300  CONTINUE 


FORM  THE  AVERAGE  CLASS  COVARIANCE  MATRIX 

SCALE  =  1.0  /  ELOAT ( NSAMP  A  NCLASS) 

DO  400  I  =  1 ,  D  IM 

DO  375  J  =  I , D I M 

DO  350  N  =  1, NCLASS 

SU(I,J>  =  SU(I,J>  +  (SCALE  A  R(I,J,N>) 

350  CONTINUE 

375  CONTINUE 

400  CONTINUE 

DO  410  I  *  1  ,DIM 

UR ITE ( DEV , A  >  (SUCIfJ),J  =  1 ,  D IM ) 

410  CONTINUE 

PUT  THE  LOWER  TRIANGLE  OF  MATRIX  SU  INTO  ARRAY  A  (FOR  SYMMETRIC 
STORAGE  MODE) 

K  =  1 

DO  450  I  *  1 , DIM 
DO  425  J  =  1 ,  I 

A ( K  >  =  SU ( I f  J ) 

K  =  K  +  1 

425  CONTINUE 

450  CONTINUE 

CALCULATE  THE  EIGENVALUES  AND  EIGENVECTORS  OF  SU  (IN  ARRAY  A) 

JOBN  =  2 
IZ  =  DIM 

CALL  EIGRS(A,DIM, JOBN,D,Z,  IZ,WK, IER) 

FORM  MATRIX  B  FROM  THE  EIGENVECTORS  /  SORT ( E IGEN VALUE ) 

IMAX  =  DIM  A  (DIM  -  1) 

DO  525  J  =■  l, DIM 

K  »  DIM  -  J  +  1 
UR  ITE ( DEV , A  ) 

UR ITE (DEV , A  >  'EIGENVALUE:  ',D(K> 

WRITE ( DEV , A ) 

DO  510  I  =  1 , D IM 

C  UR ITE ( DEV  ,  A  >  'EIGENVECTOR:  ',Z(IMAX+I) 

B (  I , J  >  =  Z ( IMAX+  I )  /  SORT ( D  <  K ) ) 

C  UR ITE ( DEV , A )  ' B (  I , J )  ' ,  I , J , B (  I ,  J  > 

510  CONTINUE 
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325  CONTINUE 

C 

C  CALCULATE  THE  AVERAGE  CLASS  MEAN  VECTOR 

C 

SCALE  =  1.0  /  ELOAT ( NCLASS ) 

DO  700  I  =  1,DIM 
SUM  =  0.0 

DO  673  N  =  1, NCLASS 

SUM  =  SUM  +  ME AN ( I, N> 

675  CONTINUE 

A VEMEAN (  I )  =  SUM  A  SCALE 
700  CONTINUE 

C 

C  SUBTRACT  THE  AVERAGE  CLASS  MEAN  FROM  EACH  CLASS  MEAN  VECTOR 

C 

DO  750  I  =  1 , D IM 

DO  725  N  =  1, NCLASS 

MEAN  <  I , N )  =  MEAN  C  I , N )  -  AVEMEAN(I) 

725  CONTINUE 

750  CONTINUE 

C 

C  FORM  THE  MEAN  SCATTER  MATRIX  SB 

C 

DO  SOO  K  =  1, NCLASS 
DO  775  I  =  1 , D  IM 

DO  760  J  =  1  ,  D  IM 

SB  <  I ,  J  )  =  SB ( I ,  J )  +  SCALE  *  (MEAN ( I ,K )  *  MEAN (  3  r  K )  ) 
760  CONTINUE 

775  CONTINUE 

800  CONTINUE 

C 

C  CALCULATE  SBSUB  =  MATRIX  SB  *  MATRIX  B 

C 

DO  900  I  =  1 , D IM 

DO  875  J  =  I, DIM 
SUM  =  0.0 
DO  850  K  =  I, DIM 

SUM  =  SUM  +  SB ( I,K )  A  B(K,J> 

850  CONTINUE 

SBSUB <  I ,  J  >  =  SUM 
875  CONTINUE 

900  CONTINUE 

C 

C  CALCULATE  SBPRIME  =  MATRIX  B  (TRANSPOSED)  A  MATRIX  SBSUB 

C 

DO  1000  I  =  1 ,  D  IM 

DO  975  J  *  1 , D  IM 
SUM  =  0.0 
DO  950  K  =  1 , D  IM 

SUM  S  SUM  ♦  B ( K , I )  A  SbSUB ( K , J ) 

950  CONTINUE 

SBPR IME ( I , J)  =  SUM 
973  CONTINUE 

1000  CONTINUE 
C 

C  PUT  THE  LOWER  TRIANGLE  OF  MATRIX  SBPRIME  INTO  ARRAY  A  (FOR 


\»wv:‘wS!* 


C  SYMMETRIC  STORAGE  MODE) 

C 

K  =  1 

DO  1050  I  =  1 , D IM 
DO  1025  J  =  1,  I 

A ( K )  =  SBPR IME (  I.J) 

K  =  K  ♦  1 

1025  CONTINUE 

1050  CONTINUE 

C 

C  CALCULATE  THE  EIGENVALUES  AND  EIGENVECTORS  OF  SBPR  IME  (IN  A) 

C 

CALL  EIGRSCA.DIM, JOBN.D.Z, IZ.UK, IER) 

C  IMAX  =  •  D  IM  A  (DIM  -  1) 

C  DO  1100  J  =  1, DIM 

C  K  -  DIM  -  J  +  1 

C  URITE(DEV.A) 

C  UR  ITE ( DEV , A )  'EIGENVALUE:  ',D(K> 

C  UR  ITE ( DEV  ,  A  ) 

C  DO  1075  I  =  1  ,  D  IM 

C  UR ITE ( DEV , A )  'EIGENVECTOR:  ',Z(IMAX  +  I) 

C 1 075  CONTINUE 

C  IMAX  =  IMAX  -  DIM 

C1100  CONTINUE 

C 

C  FORM  THE  TRANSFORMATION  MATRIX  BY  MULTIPLYING  THE  B  MATRIX 

C  TIMES  THE  FIRST  EIGENVECTOR 

C 

C  UR ITE (DEV , A  )  'THE  TRANSFORMATION  MATRIX: ' 

DO  1150  I  =  1 ,D1M 

DO  1125  N  =  1 , D IM 

DO  1115  J  =  1 ,  D  IM 

TRANS <  I ,  N  )  =  TRANS (  I , N )  +  B(I,J>  A  Z  (  J  + D IM A ( D IM-N > > 
1115  CONTINUE 

1125  CONTINUE 

C  UR  ITE ( DEV , A )  ' TRANS ( 1 ,  N ) '  , C TR ANS ( I , N ) , N= 1 , D IM > 

1150  CONTINUE 

C 

C  TRANSFORM  THE  QKIGNAL  SAMPLES  OF  EACH  CLASS 

C 

DO  1200  K  =  l.NCLASS 

DO  1175  J  =  1 , NS AMP 
DO  1160  N  =  l, DIM 
Y(J,N,K>  =  0.0 
DO  1155  I  =  1 , D I M 

Y(J.N.K)  =  Y(J,N,K)  +  ( TRANS ( I , N )  A  X C  I , J , K > > 
1155  CONTINUE 

1160  CONTINUE 

1175  CONTINUE 

1200  CONTINUE 

DELIA  =  SECNDS ( T 1 ) 

UR ITECDEV , A )  ' PREKL  EXECUTION  TIME  IS:  '.DELTA 
C  URITE(DEV.A) 'THE  TRANSFORMED  SAMPLES:' 

C  UR  ITE (DEV, 1250)  ( ( ( Y ( J, K , N ) ,K= 1 , D IM ) , J  =  1 , NSAMP ) , N= 1 , NCLASS ) 

C 1250  FORMAK 1 X , < TRD IM> F 15 . 6 ) 


RETURN 

END 


CAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA******************************** 

c 

FILE  CLASSIFY. FOR  LAST  REVISION 

AUTHOR  RUTH  D. 


JULY  24,1085 
KENNETT 


THIS  FILE  CONTAINS: 


PROGRAM  CLASSIFY 
CLSHOLE 
CLSEDGE 
CLSQB J 
OUTCONFUS 


TO  LINK: 


LINK  CLASS  1FY , COLLECT , WORK  1 , RADI , PAR AM , REG  IS , S'  SSL IBRARY : V  IDEOL IB 


C 
C 
C 
C 

c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 

CAAAAAAAAAAAAAAAAAAAAAAAAAaAAAAAAAA  A'A  aaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaa 
C  PROGRAM  CLASSIFY 

caaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaa 

c 

C  THIS  PROGRAM  CLASSIFYS  THE  OBJECTS 

c 

INTEGER  DEV, IN, OUT, STEP, INRAD, M 
INTEGER  DIN.TRDIH,  NS  A  IIP  ,  NCLASS 
INTEGER  TRN.TST.ODAT 

INTEGER  OBJECT, 0BJCT,0BJDAT(0:10, 20), NUMHOLE 
INTEGER  BOUND(10,IO),BNDINT 
INTEGER  XPIX<2500»10) ,YPIX(2500,10) 

INTEGER  CLASS<20> , EDGECLASS ( 20 > 

INTEGER  BNDCT, EDGE, OBJ, PLACE 
INTEGER  DECIS,C0NFUS<20,20> 

INTEGER  G00D<20) ,BAD<20) .COUNT 
INTEGER  TOTAL, GRDTOTAL , TOTGOOD 
REAL  QRIENT<10>,LIST(2,10,10> 

REAL  CLDIS(IO) 

REAL  TRANSC 11 ,11 ) , Y< 11 ,20) 

REAL  TA( 11 ,20) ,1(25,11,20) 

REAL  RAD (-2: 250, 10) , THETA ( 1 1 , 1 0 > , ALL ( 10 , 1 1 , 

REAL  DIS.CLSDIS, HINDIS, MAXDIS 
REAL  PERCENT 
CHARACTER  INST, SHAPE 

COMMON  DEV, IN, OUT, STEP, INRAD, M 

DATA  MAXDIS/1.7E37/ 

DATA  IN, OUT/5, 5/ 

DATA  CONFUS/ 400*0/ 

DATA  GOOD, BAD/20A0, 20*0/ 


10) 


TRN  *  l 
TST  =  2 
ODAT  *  3 
BNDINT  = 
OUTEDG  * 


OPEN ( UN IT=TRN ,  FILE* 'TRAIN.DAT'  ,  STATUS* ' OLD ' ) 
OPEN ( UN1T*TST ,  F  1LE* ' TEST . DAT ’ , STATUS* ' NEW  '  ) 
OPEN < UN IT«OD AT, FILE* 'OBJECT.DAT ' , STATUS* 'OLD ' ) 
CALL  SETUP 
CALL  THRESH 
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UR ITE ( OUT ,10) 

10  FORMATUX, ' INPUT  DIMENSION  OF  TRANSFORMED  SAMPLES:  ',*> 

READ!  IN, 15)  TRDIH 

15  FORMAT (  II ) 

UR  ITE  <  OUT , 20 ) 

20  FORMATUX,  '  INPUT  THE  SHAPE  CATEGORIC:  '  ,*) 

READ ( IN , 25 )  SHAPE 
25  FORMAT  <  A1 ) 

C 

C  READ  IN  THE  OBJECT  DATA 

C 

REAIUODAT.A)  OBJCT 
DO  100  N  =  l.OBJCT 

READ<ODAT,A>QBJ,NUMHOLE, ( OB JD AT ( K , 0 B J ) ,K=1 , NUMHOLE+ 1  ) 

OBJBAT ( 0 , OB J  )  =  NUMHOLE 
100  CONTINUE 

C 

C  READ  IN  THE  DIMENSION  OF  THE  STORED  DATA  AND  CHECK  FOR  UAL  ID  ITT 

C 

READ ( TRN  ,  A )  DIM 
M  =  DIM  -  1 

C  READ  IN  THE  STORED  DATA:  THE  TRANSFORMATION  MATRIX  AND 

C  THE  TRANSFORMED  SAMPLES  OE  THE  CLASSES 

C 

DO  110  I  =  1 ,DIM 

READ ( TRN , * )  < TR ANS C I , J ) , J= 1 , D IM ) 

110  CONTINUE 

READ <  TRN  ,  A  >  NSAMP,NCLASS 

DO  160  K  =  1 , NCLASS 

READ (TRN ,  A )  ( TA ( J , K ) , J= 1 , D IM ) 

160  CONTINUE 

URITE<TST,A)  '  TEST  SET  DATA:  ' 

URITECTST, A) 

UR  ITE ( 1ST  ,  162 )  ( 1, 1=1, DIM-1) 

FORMAK IX, 'OBJECT  HOLES  ROTATE  (  DEG )  MAXRAD',3X, 

■DIM-2) ('M',I1,7X),'M',I1,3X,' ALPHA/SQRT ( BETA > '//  ) 

'M' , II ,3X, ' ALPHA/SQRT ( BETA)  '  //  ) 

UR  ITE ( DEV  ,  A  >  '  CLASSIFICATION  DATA:' 

UR ITE<  DEV  ,  A ) 

UR ITE  <  DEV , 1 70 )  M , TRD IM , SH Ap E 
170  FORMATUX, 'AR  MODEL  ORDER:  ',12,/'  DIMENSION  OF  TRANSFORMED  ', 

A  'SAMPLES:  ',13,/'  SHAPES  CATEGORY:  ',Al,//> 

UR  ITE <  DEV, 175  >  ( J, J=1 , TRD  IM ) ,  C J , J  =  1 , NCLASS ) 

175  FORMAK IX, 'OBJECT  EDGE  '  ,  3X  ,  •  TRD IM> ( ' T ' , 12 , 7X ) , 2X , 

A  < NCLASS/ < 'B' , I 1 , 8X ) ,  'EDGE  OBJECT'/) 

COUNT  =  0 
200  CALL  VIEU 

UR  ITE ( OUT , 210)  COUNT 


210  FORMATUX,  13, '  CLASSIFICATIONS,  DONE  YET?  (TYPE  Y  FOR  YES):  ',*) 

READ (  IN , 220 )  INST 
220  FORMAT  <  A1 ) 

IF ( INST  .EQ.  'Y')  GOTO  900 

UR ITE ( OUT , 230 ) 

230  FORMATUX,  '  INPUT  THE  OBJECT  NUMBER:  ',*> 

READ  (  IN, 240)  OBJECT 
240  FORMAT  (  12 ) 

C 

C  SCAN  THE  FRAME  BUFFER  FOR  THE  EDGES  OF  THE  UNCLASSIFIED  OBJECT 

C 

T1  =  SECNDS (0,00) 

DO  250  I  =  1,10 

class:  i)  =  o.o 

250  CONTINUE 

T 1  =  SECNDS (0.00! 

CALL  SNAP 

CALL  SCAN  <  BNDCT , BOUND , XP IX , YP IX ) 

IF ( BNDC  i  .EU.  0)  GOTO  999 
CALL  INSIDE(BNDCT, BOUND) 

C 

C  CLASSIFY  USING  HOLE  INFORMATION  (IF  POSSIBLE) 

C 

CALL  CLSHOLE( OBJECT, BNDCT, BOUND, CLASS, OBJDAT, OB JCT) 

C 

C  u 1ST  THE  INNER  EDGES  IN  ORDER  OF  THE  DISTANCE  BETUEEN  THE  INNER  EDGE 

C  CENTROID  AND  THE  OUTER  EDGE  CENTROID  (LARGEST  DISTANCE  IS  FIRST  IN 

C  LIST) 

C 

CALL  ORDER(BNDCT, BOUND, LIST) 

C 

C  CLASSIFY  THE  REST  OF  THE  OBJECT  USING  THE  AUTOREGRES IVE  PARAMETERS 

C  OE  THE  OBJECT  EDGES 

C 

DO  300  OBJ  =  1, BNDCT 

IF  (  CLASS  (  OB  J  )  .  EQ  .  0)  GOTO  400 
C  UR  ITE ( 1ST ,300 )  OB JECT , BOUND ( 8 , OB J ) , OR  IENT < OB J  )  , 

C  A  INT(RAD(-1,0BJ)  ) 

C  300  FORMATUX, 2<  13, 3X)  ,F8.2,3X,  15) 

GOTO  700 

400  PLACE  =  1 

500  IF(LIST(BNDINT, PLACE, OBJ)  -NE.  0.0)  GOTO  550 

IF ( PLACE  .GT.  1)  GOTO  600 
GOTO  800 

550  EDGE  =  INT(LIST(BNDINI, PLACE, OBJ)  ) 

C 

C  CALCULATE  THE  RADIUS  VECTOR  LENGTHS 

C 

560  CALL  GETRAO(XPIX( 1 , EDGE) ,YPIX( 1 .EDGE) , EDGE, BOUND, RAD) 

C 

C  CALCULATE  THE  AUTOREGRESSIVE  PARAMETERS 

C 

CALL  PARAM(EDGE, RAD, THETA, ALL) 

UR  ITE (1ST, 575)  OB JECT , BOUND < 8 , EDGE ) , UR  IENT < EDGE ) , 

A  INT<RAD< -l .EDGE) ), (THETA < I, EDGE) ,  I* l, DIN) 

FORMAT! IX, 2 (  I3,3X>,FB.2,JX,I5,3X,  DIM-I  F9.5,F11.5> 
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nnonoo  noo 


C  CLASSIFY  THE  EDGE  USING  THE  TRANSFORMED  PARAMETERS 

C 

CALL  CLSEDGE < OBJECT, EDGE, THETA, TRANS , TA , TRD IM , NCLASS , DEC  IS) 
EDGECLASS( PLACE)  =  DECIS 
PLACE  =  PLACE  ♦  1 
GOTO  500 


CLASSIFY  THE  OBJECT  BASED  ON  ALL  THE  EDGE  CLASSIFICATIONS 


600  CALL  CLSOBJ(QBJECT, OBJ, EDGECLASS,OBJDAT,OBJCT, CLASS, 

*  B0UND<8, OBJ), TRDIM, NCLASS) 


700  DELTA  =  SECNDS(Tl) 

WRITE (DEV, 725)  OB JECT , CLASS ( OB J > , DELTA 
725  FORMAT< IX, 'OBJECT  ',12,'  IS:  ',12,/'  CLASSIFY  TIME  IS:  ', 

A  F10 . 6/ ) 


CALCULATE  THE  ORIENTATION  OF  THE  EDGES 


CALL  RQTATE(EDGE, BOUND, ORIENT) 

CALL  OUTPUT (EDGE, BOUND, RAD, ORIENT) 


750 

800 


COUNI  *  COUNT  +  1 
CONFUS ( OBJECT ,CLASS(OBJ) >  - 
IF ( CLASS ( OB J )  .EG.  OBJECT) 
BAD(OBJECT)  =  BAD ( OBJECT  >  + 
GOTO  800 

GOOD (OBJECT)  =  GOOD(OBJECT) 
CONTINUE 


CONFUS(OBJECT,CLASS(OBJ) )  +  1 
GOTO  750 
1 

♦  1 


GOTO  200 


800  CALL  OUTCONFUS ( CONFUS , OB JCT ) 


UR ITE ( DEV , A  ) 

UR ITE ( DEV  ,  A  )  'CLASS  PERCENT  CORRECT' 

GRDTOTAL  =  0.0 
DO  850  K  =  1, OB JCT 

TOTAL  =  GOUD(K)  ♦  BAD ( K ) 

IF ( TOTAL  .EO.  0)  GOTO  950 

PERCENT  =  100.0  A  ELOAT ( GOOD ( K ) )  /  FLOAT ( TOTAL  > 
UR ITE<  DEV , 940 )  K, PERCENT 
940  FURMAK IX, I3,5X,F10.5> 

GRDTOTAL  =  GRDTOTAL  +  PERCENT 
950  CONTINUE 


PERCENT  *  GRDTOTAL  /  FLOAT ( OBJCT ) 

UR ITE( DEV ,960 )  PERCENT 

960  FORMAT< IX, 'OVERALL  PERCENT  CORRECT:  ',F10.5> 
GOTO  1000 


999  UR ITE < OUT, A )  'NO  OBJECTS  FOUND' 

1000  CLOSE(UNIT-ODAT) 

CLOSE(UNIT«IRN) 

CLOSE ( UN IT»TST > 

STOP 

END 


non  ono  ooooooo 
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CLSHOLE 

SUBROUTINE  CLSHOLE ( 08 JECT , BNDCT , BOUND , CLASS , 08 JD AT , OB JCT > 

THIS  SUBROUTINE  CLASSIFIES  THE  OBJECT  ON  THE  BASIS  OF  THE  NUMBER 
OF  HOLES  IT  HAS. 


INTEGER  DEV, IN, OUT 

INTEGER  BNDCT, BOUND <10, 10 ), CLASS ( 10) 

INTEGER  BND , OBJ 

INTEGER  MAXHOLE , NUMHOLE, NHOLE (2,10) 

INTEGER  OB J ECT, OB JCT, OB JDAT< 0:10,20) 

COMMON  DEV,  IN, OUT 

DATA  NH0LE/20 AO/ 

BND  =  1 
OBJ  =  2 
MAXHOLE  =  10 

00  100  I  »  1,10 
DO  50  J  =  1,2 

NHOLE ( J , I )  =  0 
50  CONTINUE 

100  CONTINUE 

DO  400  N  »  1, BNDCT 

DO  300  K  =  1 ,  OBJCT 

IF  <  BOUND  <  8 ,N)  .NE.  0BJDAT<0,K)>  GOTO  300 
NUMHOLE  »  BOUND  <  8 , N ) 

IF(NHOLE<BND, NUMHOLE)  .EQ.  0)  GOTO  250 
NHOLE<BND, NUMHOLE)  -  -1 
GOTO  300 

250  NHOLE<BND, NUMHOLE)  =  N 

NHOLE (OBJ, NUMHOLE)  =  K 
300  CONTINUE 

400  CONTINUE 


DO  500  NUMHOLE  *  l, MAXHOLE 

IF(NHOLE< BND, NUMHOLE)  .LE.  0)  GOTO  500 
CLASS(NHOLE(BND, NUMHOLE) )  «  NHOLE ( OB J , NUMHOLE ) 

UR  HE (DEV, 450)  OB JECT , NHOLE <  BND , NUMHOLE ) , NUMHOLE, 

A  NHOLE(OBJ, NUMHOLE) 

50  FORMAT* IX, 2< 13, 3X> , 'CLASSIFICATION  BASED  ON:  ',11,'  HOLES:', 13) 

500  CONTINUE 

CLOSE(UNII-LUN) 

RETURN 
END 


CLSEDGE 


SUBROUTINE  CLSEDGE < OB JECT , EDGE , THETA , TRANS , TA , TRD IM , NCLASS , DEC  IS ) 


non  nnn  n  o  n  o  o  o 
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THIS  SUBROUTINE  CLASSIFIES  AN  EDGE  OF  THE  OBJECT.  THE  CLASSIFICATION 
IS  BASED  ON  THE  SMALLEST  DISTANCE  FROM  THE  TRANSFORMED  PARAMETERS  OF 
THE  UNCLASSIFIED  EDGE  TO  THE  AVERAGE  OF  THE  TRANSFORMED  SAMPLES  OF  THE 
TRAINING  SET  OF  EDGES. 

INTEGER  DEV, IN, OUT, STEP, INRAD, M 
INTEGER  EDGE,TRDIM,NCLASS,DECIS 


INTEGER 

OBJECT, MINCLS 

REAL 

MAXDIS, HINDIS 

REAL 

THtTA(ll,10>,TA(ll,20) 

REAL 

TRANS (11,11), Y(ll, 20) 

REAL 

CLSDISf 10) 

COMMON 

DEV,  IN, OUT, STEP, INRAD  ,  M 

DATA 

MAXDIS/1.7E37/ 

N  =  EDGE 
DIM  =  M  ♦  1 

TRANSFORM  THE  AR  PARAMETERS  OF  THE  UNCLASSIFIED  EDGE 

DO  200  J  =  l.TRDIM 
Y(J,N>  =  0.0 
DO  ISO  I  -  1 , D IM 

Y  ( J ,  N )  *  Y  ( J ,  N  >  ♦  ( TRANS ( 1,3)  A  THETA ( I , N ) ) 

ISO  CONTINUE 

200  CONTINUE 

FIND  THE  SMALLESI  DISTANCE  TO  A  CLASS 

HINDIS  *  MAXDIS 

MINCLS  =  0 

DO  400  K  *  I , NCLASS 
CLSDIS(K)  =  0.0 
DO  300  J  =  l.TRDIM 

CLSDIS(K)  =  CLSDIS(K)  +  (Y(J,N>  -  TA(J,K))AA2 
300  CONTINUE 

IF(CLSDIS(K)  .GT.  M1NBIS)  GOTO  400 
HINDIS  =  CLSDIS(K) 

MINCLS  =  K 

400  CONTINUE 

DEC  IS  =  MINCLS 

WRITE (DEV, 450)  OB JECT , EDGE , ( Y ( J , N ) , J  =  1 ,  TRD  IM  )  , ( CLSD IS ( K ) , K= 1 , NCLASS  )  , 
A  DECIS 

450  FORMAT < 1 X , 2  < 13 , 3X ) , •  TRD IM>F 1 2 . 7 , •  NCLASS; F 1 2 . 7 , 8X , 13 > 

RETURN 

END 

C - 

C  CLSOBJ 

C - - - - - 

SUBROUT INE  CLSOB J ( OB JECT , OB J , EDGECLASS , OB JDAT , OBJCT , OB JCLASS , THOLE , 

A  TRDIM, NCLASS) 

C 


nnnn  o  o  o 
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C  THIS  SUBROUTINE  CLASSIFIES  THE  OBJECT  USING  THE  EDGE  CLASSIFICATIONS 


INTEGER  DEV, IN, OUT 

INTEGER  EDGECLASSf 10) ,0BJDAT(0: 10,20) 

INTEGER  OBJECT, OBJ, OBJCT,OBJCLASS(10) 

INTEGER  COUNT, OBJHOLE, THOLE 
INTEGER  TRD  IN , NCLASS 

COMMON  DEV, IN, OUT 

N  =  1 

100  IF  <  N  .GT.  OBJCT)  GOTO  600 

OBJHOLE  =  OB JDAT ( 0 , N  > 

IF ( THOLE  .N£.  OBJHOLE)  GOTO  300 

COUNT  =  0 

DO  200  J  =  1 , OB JHOLE+ 1 

IF ( EDGECLASS ( J  >  .NE.  OBJDAT<J,N>>  GOTO  200 
COUNT  =  COUNT  +  1 
200  CONTINUE 

IF ( COUNT  .Ed.  OBJHOLE+1)  GOTO  500 

300  N  =  N  +  1 

GOTO  100 

500  OB JCLASS ( OBJ )  =  N 

GOTO  1000 

600  OB JCLASS ( OB J  )  =  0 

1000  CONTINUE 

C  1000  WRITE (DEV, 1010)  QB JCLASS ( OB J ) 

C  1010  F0RMAT(1X,12X,<TRD IM> (8X ),<NCLASS>(10X),14X,I3) 
RETURN 
END 


OUTCONFUS 


SUBROUTINE  OUTCONFUS < CONFUS , NCLASS ) 

THIS  SUBRUTINE  WRITES  OUT  THE  CONFUSION  MATRIX  OF  THE  CLASSIFI¬ 
CATION  DATA 

INTEGER  C0NFUS(20, 20) .NCLASS 
INTEGER  DEV 

COMMON  DEV 

UR ITE ( DEV , 100 ) 

100  FORMAT ( IX , ' THE  CONFUSION  MATRIX  OF  THE  CLASSIFICATION  DATA'//) 

UR  ITE (DEV , 150 )  (1,1*1 .NCLASS ) 

150  FORMAT ( 1X.6X, <NCLASS> II»,/7X, (NCLASS A 5> ( '-' ) ) 

DO  400  I  *  1, NCLASS 

UR  ITE (DEV, 200)  I, (CONFUS ( I , J >, J= 1 , NCLASS > 

200  FORMAT( IX, 15, ' | ' ,<NCLASS> 15) 

400  CONTINUE 

tCTU«H 

w* 
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CAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA 
FILE  COLLECT. FOR  LAST  REVISION  JUNE  10,  1985 


SUBROUTINES  CONTAINED  IN  THIS  FILE  ARE: 


INSIDE 

ROTATE 

OUTPUT 
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INSIDE 


SUBROUTINE  INS  IDE  (  BNDCT , BOUND  ) 

THIS  SUBROUTINE  FINDS  OUT  UHICH  BOUNDARIES  ARE  INSIDE  OTHER  BOUNDARIES 

INTEGER  BNDCT , BOUND! 10,10) , BND INT 
INTEGER  XR,XL,XPREV,RPIXEL,Z 
INTEGER  DEV, IN, OUT, STEP, INRAD 
INTEGER  OBJ, COUNT 

INTEGER  YLINE, FIRSTX, LASTX, STARTX, XINC,  BUFF ( 0:511) 

LOGICAL  EDGE,RBND(20),LBND(20) 

COMMON  DEV, IN, OUT, STEP, INRAD 

DO  40  I  -  1, BNDCT 
BOUND (7,1)  »  0 
BOUND ( 8 , I)  -  0 
40  CONTINUE 

FIRSTX  =  0 
LASTX  =*  511 
XINC  *  l 
EDGE  ”  .TRUE. 

LOOP  FOR  EACH  EDGE 


DO  600  N  =  1, BNDCT 

BND INT  =■  BOUND ( I , N  > 

STARIX  =  BOUND ( 2 , N ) 

YLINE  =  BOUND ( 3 , N  > 

DO  100  I  =  1, BNDCT 

RBND(l)  =  .FALSE. 

LBND  < I )  =  .FALSE. 

100  CONTINUE 

SEARCH  FOR  EDGES  ON  THE  RIGHT 

XPREV  =  0 

CALL  RSC AN (YLINE, FIRSTX, LASTX, XINC, BUFF) 

DO  200  XR  «  STARTX, LASTX 
Z  *  BUFF  <  XR ) 

IF  < (Z.EQ.O) .OR. (Z.E0.255) .OR. ( Z . ED . BND INT ) )  GOTO  200 


noon  non  non  non 


<  XPREV  ♦  1))  GOTO  150 


IF ( XR  .NE. 

XPREV  »  XR 
GOTO  200 

150  RBNO(Z)  »  RBND  (  Z  >  .NEQV .  EDGE 

XPREV  =  XR 
200  CONTINUE 

SEARCH  FOR  EDGES  ON  THE  LEFT 

XPREV  =  0 

DO  300  XL  =  STARTX,FIRSTX,-1 
Z  =  BUFF ( XL ) 

IF  <  (Z.EQ.O)  .OR.  (Z.EQ.255)  .OR.  (  Z  .  EQ .  BND  INI )  )  GOTO  300 
IF ! XL  .NE.  (XPREV  -  1))  GOTO  250 
XPREV  *  XL 
GOTO  300 

250  LBND(Z)  *  LBND(Z)  .NEQV.  EDGE 

XPREV  =  XL 
300  CONTINUE 

COMPARE  THE  RIGHT  AND  LEFI  TALLIES 

COUNT  =  0 

DO  400  I  =  1, BNDCT 

IF < . NOT . <  RBND ( I >  .AND.  LBND(I)>)  GOTO  400 
COUNT  »  COUNT  ♦  1 
OBJ  «  I 

BOUND  (  7 ,  N )  »  1  ! FOUND  AN  OUTER  EDGE 

400  CONTINUE 

IE (COUNT  .EO.  1>  GOTO  450 
BOUND ( 7 , N )  »  0 
GOTO  GOO 

450  B0UND(8,0BJ)  =  B0UND(8,0BJ)  +  1  !COUNT  HOLES  INSIDE  OUTER  EDGE 

600  CONTINUE 

RETURN 

END 


ROTATE 


SUBROUTINE  ROTATE ( BNDCT , BOUND , OR IENT > 

THIS  SUBROUTINE  CALCULATES  THE  ROTATION  OF  THE  MAXIMUM  RADIUS  VECTOR 
(IN  DEGREES)  IN  REFERENCE  TO  THE  STANDARD  CARTES  ION  COORDINATE  SYSTEM. 

INTEGER  BNDCT, BOUNDdO, 10) 

INTEGER  XC,YC,XMAX,YMAX 
REAL  OR  IENT ( 10 ) 

REAL  XDIF,YDIF, THETA 

XC  -  4 

YC  »  5 
XMAX  *  9 
YMAX  »  10 

DO  500  I  »  1, BNDCT 

XD IF  «  FLOAT ( BOUND < XMAX ,  I)  -  BOUND(XC,I>> 

YD  IF  *  -FLOAT( BOUND* YMAX,  I)  -  80UND(YC,I)> 
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IE ( XD IE  .NE.  0.0>  GOTO  400 
IE<  YD  IE  .GT.  0.0)  GOTO  300 
THETA  =  90.0 
GOTO  490 

300  THETA  =  270.0 

GOTO  490 

400  THETA  =  AIAND(YDIF/XDIF> 

IE ( XD IE  .LI.  0.0)  THETA  =  THETA  h 
IF ( THETA  .LI.  0.0)  THETA  =  THETA 
490  ORIENT(I)  =  THETA 

500  CONTINUE 

RETURN 

END 
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SUBROUTINE  OUTPUT (BNDCT, BOUND, RAD, ORIENT) 

THIS  SUBROUTINE  PRODUCES  THE  OUTPUT  OF  THE  OBJECT  DATA 

INTEGER  BNDCT, B0UND<10, 10) 

INTEGER  DEO 

REAL  RAD(-2:300,10),0RIENT(10) 

COMMON  DEV 
WRITE (DEV, SO) 

FORMAT (IX, 'EDGE' ,6X,' START' ,8X, 'CENTER' ,6X, 'NUMPIX' ,2X, 

'  INSIDE' ,2X, '  HOLES','  ROI  (DEG)','  #  RAD  VEC', 

'  MAX  RAD  VEC  (PIX) '/) 

DO  800  N  =  1, BNDCT 

WRITE (DEV, 750)  ( BOUND (  I,N> , 1  =  1 ,8) , OR 1ENT < N ) , RAD ( -2 , N ) , 
RAD ( -  1 ,N  > 

FORMAT ( IX, 13, 5X, ' ( ' , 13, ' , ' , 13, ' ) ' ,5X, ' ( ' , 13, ' , ' , 

13, '  )  ' ,4X,  I4,2(5X, 13 > ,GX,FS. 1 ,2( 4X,FS.2> ) 

CONTINUE 

WRITE(DEV,810) <N,N  =  1, BNDCT) 

FORM AT ( / / , <BNDCT>  (  6X ,  'EDGE'  >/<BNDCTX7X,  '(',  II ,')'),/ ) 

UR ITE ( DEV , 320 ) 

FORMAK/,'  NUMBER  OF  RADIUS  VECTORS:',/) 

UR  ITE (DEV, 825)  (RAD(-2,N> ,N=1, BNDCT) 

F0RMAT(1X,3X,8E10.4> 

UR  ITE ( DEV , 330 ) 

FORMAK/,'  MAXIMUM  RADIUS  VECTOR:  ',/) 

UR  ITE (DEV, 835)  ( RAD ( -1 , N ) , N= 1 , BNDCT > 

FORMAK 1X,3X,8F10 . 4 ) 

WRITE ( DEV , 850 ) 

FORMAK/,'  RADIUS  VECTOR  LENGTHS:'/) 

DO  950  1=1, INRAD+1 

UR  ITE< DEV , 920 )  I, ( RAD( I , N> , N= 1 , BNDCT ) 

FORMAT ( IX , 13 , 8F10 . 4  > 

CONTINUE 


no  no 


3X0  IF(LIST(DIS,PLACE+1 ,0BJ>  .EQ.  0.0)  GOTO  320 

PLACE  =  PLACE  *■  1 
GOTO  310 

320  LIST(BNDINT,PLACE+1,0BJ>  =  L 1ST < BND INI , PLACE, 0B J ) 

LISI<DIS,PLACE+1,0BJ>  =  LISX(DIS, PLACE, OBJ) 

IE (PLACE  .LE.  HARK)  GOTO  340 
PLACE  =  PLACE  -  1 
GOTO  320 

340  LISKBNDINT, PLACE, OBJ)  =  FLOAT  ( N  > 

LIST<DIS, PLACE, OBJ)  =  DIST 

500  CONTINUE 

UR  ITE  <  DEO , 550 ) 

550  FORMAT( IX,//'  OUTER  EDGE  INNER  EDGE  ORDER  DISTANCE'  /> 


720 


730 


740 


800 


DO  800  OBJ  =  1 , BNDCT 
PLACE  =  1 

IF(LIST(BNDINT, PLACE, OBJ)  .EQ.  0.0)  GOTO  800 
INEDGE  =  INT(LIST(BNDINT, PLACE,  DBJ) ) 

IF(  INEDGE  .EQ.  OBJ)  GOTO  730 

QRD  »  PLACE 

GOTO  735 

INEDGE  =»  0 

ORD  ■  0 

CONT INUE 

UR  ITE ( DEU , 740 )  OB J , INEDGE , ORD , L 1ST ( D IS , PLACE , OB J > 
FORMAT (IX, 16, 9X, 13,1 IX, I2,7X,F10.4> 

PLACE  *  PLACE  ♦  1 
GOTO  720 
CONTINUE 


RETURN 

END 


